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ABSTRACT 


The  problem  of  identifying  the  error  parameters  associated  with  inertial 
measurement  units  is  considered  in  this  report.  This  is  an  important 
practical  problem  which  is  included  in  a  large  class  of  system  param¬ 
eter  identification  problems.  A  general  approach  for  formulating  the 
many  possible  inertial  measurement  unit  (IMU)  error  parameter  configur 
ations  is  given,  and  specific  realizations  are  specified  in  detail.  Tne 
formulation  is  such  that  time  correlated  environmental  and  observational 
random  disturbances  can  be  incorporated.  Experimental  results  showing 
the  effects  of  state  and  observation  noise  power  levels,  and  the  nominal 
trajectory  on  the  identification  of  the  error  parameters  for  three  specific 
configurations  are  presented.  These  results  indicate  that  a  meaningful 
optimization  problem  can  be  formulated  in  terms  of  the  nominal  trajectory 
variables.  The  problem  is  then  considered  as  an  optimal  control  problem 
with  the  cost  being  a  functional  of  the  estimation  covariance  matrix  and  the 
controls,  where  certain  nominal  trajectory  variables  are  considered  as  the 
controls.  The  question  of  the  existence  of  optimal  controls,  the  necessary 
conditions  which  the  optimal  controls  must  satisfy,  and  the  computational 
aspects  for  computing  the  optimal  controls  are  considered. 
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NOMENCLATURE 


En  n-dimensional  Euclidean  space 
x,x  n-dimensional  vectors 

A, {A1  nxm-dimensional  matrices 

x*,A*  transpose  of  a  vector  and  a  matrix,  respectively 
2 

<x,y>,  f|x!|o  inner  product  of  two  vectors 

t,  t0,T,I  time,  initial  time,  final  time,  and  the  time  interval 
[to,T],  respectively 

d>(t,r)  fundamental  solution  (transition  matrix)  ol  a  linear 
differential  equation 

IMU^  j  inertial  measurement  unit  (•) 

object  j  in  the  system,  with  j  =  m.o,s,r,  and  i, 
respectively,  to  denote  the  vehicle,  a  reference 
IMU  (IMUq)  ,  a  slaved  IMU  (IMUg),  a,  generally, 
non-stationary  reference  point,  and  an  inertialiy 
fixed  reference  point,  respectively 

B.  P.  ^  a  fixed  point  of  object  j 

,  <*-  1, 2, 3,  orthogonal  axes  centered  at  B.  P.  ^ 

b(J»k),  €b<i  »k)  vector  distance  from  B.P.  O)  to  B.P.  and 

the  error  in  the  knowledge  of  b(j>k),  respectively 

a,  [&>]  angular  velocity  vector  and  matrix  respectively 

a,  [a]  acceleration  vector  and  matrix,  respectively 

[nun1]  especially  defined  angular  velocity  matrices 

e0»k)  a=  i,2,3,  body  angles  from  to  X^k) 

k)  orthogonal  coordinate  transformation  from  B^  to  B*k) 


n  random  position  changes 


x  a  time  derivative  with  respect  to  a  non-stationary  coordinate 
system 

Aa^^  difference  in  sensed  acceleration  between  and  3^) 

ABm,k^  difference  in  measured  acceleration  between  andB^) 

vjtS’ a\  Pm’ first  second  integrals,  respectively,  of  the 
measured  acceleration  between  IMU0  and  IMtL 

o  O 

n,n,  v,_pc,^w  various  random  processes 

*>*(*0)'  *v  i  =  1>2>3  vector  of  initial  misalignment  angles 

Ka  vector  of  the  (six)  accelerometer  error  parameters 

Kg  vector  of  the  (six)  gyroscope  error  parameters 

Km  vector  of  various  measurement  error  parameters 

~S»  kS->  3,  mass  unbalance  errors  along  the  spin  axis  for 

1  tiiree  gyros 

kIj,  i  =  l,2 , 3,  mass  unbalance  errors  along  the  input  axis  for 
three  gyros 

£,  «j,  i  =  1,2,3,  constant  drifts  for  three  gyros 

Sjj^R^t)  power  spectrum  and  autocorrelation,  respectively, 
for  the  random  process  x 

Configuration  IA  Identification  of  £,  with  correlated  random 
acceleration 

Configuration  IB  Identification  of  ♦ ,  with  uncorrelated  random 
acceleration 

Configuration  n  Identification  of  ♦  and  £ 

Configuration  QI  Identification  of  and  the  mass -unbalance  gyro 
drUts'  *m.u. 

cjj  =^0*  U>e  correlated  random  acceleration  power 
2 

a  r*I  the  uncorrelated  random  acceleration  covariance  matrix 


9 

o-^-i  the  observation  noise  covariance  matrix. 

.  the  standard  deviation  in  the  estimated  misalignment  angles 

(j/\_  the  standard  deviation  in  the  estimated  constant  gyro 
€i  drift  rates 

the  standard  deviation  in  the  estimated  mass -unbalance 
1  1  gyro  drift  terms 

Efx],  Cov(x,  expected  value  of  x,  and  covariance  of  x  and  y, 
respectively 

P  =  E((x  -£)(x  -£)*]  covariance  matrixof  errors  in  the  estimated 

value  $  from  the  true  value  x 

p..  elements  of  the  covariance  matrix  P 
ij 

P3^  covariance  matrix  associated  with  the  x  and  y  portions  of  a 
vector  z,  where  z  =  (x,  y)* 

!l  P  Jlj ,  i  =  1, 2, 3,™,  various  norms  in  En, associated  with  matrix  P 

R,  Q  covariance  matrices  of  the  time  uncorrelated  observation 
and  state  vector  noise  disturbances. 

u,fx  an  m-vector  of  control  functions  and  a  scalar-valued  control 
function,  respectively 

Mj  various  constraints  on  the  control  function  u 

U(t)CEm  a  pointwise  constraint  set  for  u(t)  for  t«  I 

a  set  of  control  functions 

J  a  cost  functional  associated  with  P,u,  and  T 

Trace  |WP(Tiu)|  sum  of  the  diagonal  elements  of  WP  evaluated 
L  J  at  t  =  T,  with  the  control  u(t),  tQ<  tsT 

{ ujj}  a  sequence  cf  functions 

Ujj  -2-  u  weak  convergence  of{u^)  tou 

L2([a,  b) )  space  of  square  integrable  vector  funcuons  on  the 
time  interval  (a,  bj 
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H  Hamiltonian  function  associated  with  P  u  and  J 

A,  Xy,  vector  of  adjoint  variables  associated  with  the  ©variance 
matrix  P,  and  the  elements  of  A,  respectively 

I(A0),  I  (IV)  cost  functional  associated  with  the  error  in  guessing 
the  initial  adjoint  matrix  A^ tg)  =A0>  and  in  guessing 
the  final  covariance  matrix  P(T)  =PT,  respectively 

W  a  non-negative  weighting  matrix  associated  with  P. 


CHAPTER  I 


INTRODUCTION 

The  problem  of  Identifying  the  error  parameters  associated  with 
inertial  measurement  units  (IMU)  is  considered  in  detail.  This  is 
an  important  practical  problem  which  is  included  in  the  general 
class  of  problems  formulated  in  the  next  section.  A  general 
approach  for  formulating  the  many  possible  IMU  error  parameter 
configurations  is  given  in  Chapter  II,  and  specific  realizations  are 
given  in  Chapter  m.  Minimum  variance  estimation  techniques  for 
the  error  analysis  of  the  identification  procedure,  and  typical  ex¬ 
perimental  results  are  presented  in  Chapters  IV  and  V.  In  Chapter 
VI,  techniques  from  optimal  control  theory  are  used  to  characterize 
the  best  nominal  trajectory  from  a  class  of  admissible  nominal 
trajectories.  This  trajectory  is  best  in  the  sense  that  the  identi¬ 
fication  of  the  error  parameters  is  the  mcst  satisfactory  under  given 
physical  conditions.  Computational  aspects  for  obtaining  the  optimal 
nominal  trajectory  are  discussed  in  Chapter  VII. 
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1. 1  The  General  Problem 


Many  physical  systems  might  be  described  by  the  state  equation 

y  =  f(t,  y,  u,}  w, ,  a, )  (1.1) 

and  the  observation  equation 

z  =  h(t,  y,  w2,  )  +  v  (1.2) 

In  these  equations  all  the  variables  are  generally  vectors,  except 
for  the  time  t .  The  n-dlmensional  vector  function  f  and  the 
m-dimensional  vector  function  h  are  assumed  to  be  sufficiently 
smooth  so  that  the  required  linearizations  which  are  subsequently 
required  are  valid.  The  specific  definitions  are  as  follows: 

y  is  an  nxl  vector  of  state  variables, 
u,  is  an  r,  x  1  vector  of  state  control  variables, 

Wj  is  a  p,xl  vector  of  random  state  disturbances, 

arj  is  a  k,xl  vector  of  constant  state  parameters, 

z  is  an  in  xl  vector  of  observation  variables, 

u2  is  an  r2x  1  vector  of  observation  control  variables, 

w*  is  a  BgXl  vector  of  random  observation  disturbances, 

a2  is  a  kfcx  l  vector  of  constant  state  parameters, 

v  is  an  mxl  vector  of  random  white  noise  disturbances. 
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The  vectors  and  a2  of  system  parameters  are  assumed  to  be 
made  up  of  constants,  although  smooth  functions  would  also  be 
considered  in  this  framework  by  making  suitable  polynomial 
approximations,  and  then  increasing  the  dimension  of  and  a2 
appropriately.  The  random  disturbances  Wj  and  w2  may  be 
time-correlated  or  not,  but  it  is  assumed  that  the  statistics  of 
these  vectors  are  completely  known  and  that  they  may  be  modelled 
by  white  noise  passing  through  a  linear  dynamical  system 
(Markovian  property) . 

The  general  problem  which  might  be  presented  is  to  identify  the 
system  parameters  and  cv2  based  on  the  observations  z.  Of 
course,  it  is  also  usually  of  interest  to  obtain  estimates  of  the 
state  y.  The  controls  uj  and  u2  are  functions  which  are  avail¬ 
able  to  change  the  evolution  of  the  state  and  observation  equations. 
This  change  in  these  equations  is  to  be  made  in  such  a  way  that  the 
identification  might  be  accomplis^l  in  some  optimal  fashion, 
subject  to  certain  constraints  in  the  controlling  functions  u^  and  u2. 
Except  for  special  theoretical  questions,  there  is  not  much  that 
can  be  said  about  the  analytical  properties  of  this  problem  in  this 
generality.  Certain  linearizations  might  be  made  to  make  the 
problem  more  tractable  and  perhaps  ame^-ble  to  certain  mean 
square  optimization  techniques. 


Under  suitable  assumptions,  the  system  (l.l)  ,  (1.2)  may  be 
written  as 


y  =  Ay(t , Uj )  y  +  Bwl(t,Uj)  w,  +  CQ>  (t, u^a,  (1.3) 
and 

z  =  My(t,u2)y+BW2(t,u2)  w2  +  Ca2(t,u2)a2+  v 

(1.4) 

These  equations  represent  small  deviations  about  nominal  or 
expected  values  of  the  variables  described  above,  which  are 
denoted  by  y,  w,  etc.  The  matrices  in  Equations  (1.3)  and  (1.4) 
are  the  various  partial  derivatives  of  the  functions  f  and  h,  and 
are  evaluated  about  the  nominal  or  expected  values.  Typically, 


and 


Ay(t,  u,) 


af(t,  • ) 

ay 


BW,(t,u,)  H  af-fa  •> 


It  is  noticed  that  we  do  not  linearize  about  the  controls  uj  and  u £ 
because  in  the  application  which  is  considered  here,  the  nominal 
values  flj ,  CI2  are  also  to  be  chosen  in  an  optimal  way.  It  is 
next  assumed  that  ,  W£  have  the  Markovian  representations 
Wj  =  Aw*  Wi  +  BW|'  v1  and  w2  =  AWz  w2  +  Bw2  v2  * 
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If  we  write  <*j  =  0  and  =  0,  then  the  above  equations  might 


be  collected  together  and  written  as 


i  ~  AY  +  Bw 


(1.5) 


z  »  MY  +  v 


(1.6) 


where 


Y  s  (y,  w,,  w2,ot}  ,a2) 


Ay  BW]  0  C®1 


0  AW1  0 


A  =  A(t,uj)»  0  0  Aw2  0 


w  =  (w,,  w2) 


Bw1*  0 
0  BW2 


« M(t,u2)=r 


I  0  Bw2 


0  (f^J 


If  the  controls  Uj ,  U£  have  been  prespecified,  the  problem  formu¬ 
lated  by  equations  (1. 5)  and  (1.6)  is  in  the  standard  form  for 
applying  the  well-known  minimum  variance  estimation  techniques 
[l].  Assuming  the  random  vectors  possess  the  known  statistical 


properties 


(w(t)J  =  E  [w(r)J  =  0,  E  lw(t)  v(r)*)  a  0, 


Efw(t)  w(r)V=Q(t)6(t-7),  E[v(t)  v(rf]  =  R(t)6(t-r), 


E(Y(t0)J  =  0,  EiY(t0)Y(t0f]  =  P0, 
EIY(to)v(r)*!  =  0,  ErY(t0)w(T)*]  =  0, 


then  the  minimum  variance  estimate  Y  is  given  by  the  solution  of 


the  differential  equation 


=  [a  -km]$ 


+  Kz  ;  Y(t0)  =  0 


(1.7) 


The  optimal  gain  matrix  K  is  given  by 


K  =  PM*R 


where  the  covariance  matrix  P,  which  is  defined  by 
Trace  P  =  Trace  eJc?  -  Y)(Y  -  f  fj, 


satisfies  the  matrix  Riccati  equation 
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(1.8) 


dP/dt  =  AP  +  PA*  -  PM*  R‘XMP  +  BQB*  ; 

P(t0)  =  E[y(t0)y(t0)*l  =  P0 

The  dependence  of  the  solution  on  the  control  parameters  us(uj,u2) 
will  be  denoted  by  and  Pu .  The  optimal  control  problem  which 
naturally  arizes  is  to  minimize  some  functional  of  P,  for  example 
njin  tracefWPfTjuJ]  where  W  is  a  weighting  matrix  and  T  is  the  end 
time,  or  else  to  choose  a  u  such  that  certain  diagonal  elements  of 
P  fall  into  a  prespecified  region  in  minimum  time. 

For  a  linearized  minimum  variance  estimator  this  cost  would 
depend  on  a  particular  stochastic  realization  of  the  system.  The 
variances  in  the  resulting  P(T)  for  particular  realizations  would 

i 

then  need  to  be  considered  to  decide  on  a  suitable  ’'average”  cost. 

Thus  a  linearized  minimum  variance  estimator  would  complicate 
the  optimization  considerably.  Since  we  should  like  to  apply  some 
of  the  known  techniques  from  optimal  control  theory  D3.  it  shall  be 
assumed  that  we  have  a  linear  system  from  the  outcet.  This  is 
indeed  the  case  for  the  IMU  error  identification  problem,  because 
the  design  of  the  associated  instruments  is  such  that  it  keeps  all 
errors  reasonably  small.  So  that  even  though  the  various  errors 
might  be  very  complicated  in  nature,  their  contribution  to  the  system 
dynamics  can  be  accurately  represented  by  a  linear  transformation. 
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ssi'-a  fcf  - 


1.2  The  2MU  Error  Parameter  Identification  Problem 


The  problem  of  identifying  the  error  parameters  associated  with 
inertial  measurement  units  (IMITs)  arizes  in  many  applications 
associated  with  the  inertial  navigation  of  ships,  airplanes,  missiles, 
and  space  vehicles.  Investigation  of  this  specific  class  of 
•  important  problems  provides  motivation  and  insight  into  the  general 
problem  outlined  above.  A  discussion  of  the  error  parameters 
which  are  typically  associated  with  IMU’s  and  the  physical  en¬ 
vironment  in  which  such  systems  operate  is  given  in  Chapter  II. 

A  general  formulation  of  the  problem  and  a  discussion  of  experi¬ 
mental  results  is  given  in  Chapters  II  to  V 

Generally,  the  IMU  error  identification  problem  has  the  form 

*  =  Ax  x  +  Bn  n+  Ca<*+  wx 
ft  =  An  n  +  wn 

(1.9) 

6=0 
Z  =  X+  V 

Equations  (1.7)  and  (1.8)  then  become,  respectively, 

d3/dt  =  (Ax  -P^R"1)  k  +  Bn  k  +  Ca&  +  P^R’1  z 
dfi/dt  =  -P^R*1  x  +  Anfi  +  P^z  (1.10) 

dtt/dt  =  k  +  P^lf1  z 


and 


d  p*3^  Ax  j>xx + pxx  Ax*+  gn  pF+  pxnEp*+caIpx+  p^+qf* 

dt 

h  „nn  tnjin  „nn  4n  nx  -l_xn  _nn 
JLp  s  AP  +P  A  -P  R  P  +Q 
dt 

JLiT-  -p^r  1pxct 

dt  (i.u) 

_SL  pm=  Axpx»+PxnAn‘+  BnPnn+<fp“n-PxxR’1Pxn 
dt 

p*“s  AXPXa  +  BnPn“+C0fPp!a  -PXXR-1PXOf 
dt 

h  _  lior  .  n_na  _nx  -1  xa 
-*r  P  =  A  P  -P  R  P 
dt 

where  the  covariance  matrix  P  has  been  partitioned  in  the  obvious 
way 


"XX 

„xn 

xoT 

p 

P 

P 

pnx 

pnn 

pna 

ax 

P 

p“" 

Generally,  the  control  parameters  would  appear  in  the  matrices 
Ax ,  Bn ,  Ca,  and  Q30*.  Thus  it  is  difficult  to  study  the  optimization 
problem  in  detail  at  this  level  of  generality,  although  the  optimi¬ 
zation  problem  can  be  clearly  stated.  If  the  assumptions  used  for 
special  configurations  (configurations  IB,  II,  and  HI  of  Chapter  HI) 
are  invoked,  the  optimization  starts  to  be  tractable.  Briefly, 
these  assumptions  are  n  =  0  (no  correlated  process  noise),  and 
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diagonal  noise  covariance  matrices  Q**  *  •  I  and  R  =  cr^  •  I . 

Equation  (1. 10)  and  (1.11)  reduce  to 


-f-  &  =  -^pXX$  -  z]  +  Ca& 
dt  w 

£ a  -  *. 


(1.12) 


and 


jLp“.  c“r^+  p^c0*  -^p^p301  +  ^  •  i 
£*»“-  (1.13) 

fpF=  C«P““ 
dt  w 


respectively.  The  trajectory  control  parameters  are  now  associa¬ 
ted  with  the  matrix  Ca  only.  The  optimization  problem  is  then  to 
choose  these  trajectory  control  variables  so  that 

Trace  [w“p”(T)j, 

where  W°  is  a  positive  definite  weighting  matrix,  is  minimized  at 
some  gben  end  time  T,  or  else  to  choose  the  control  variables  such 
that 

pjj(T)<  ,  «j  specified,  |-i,***,  A; 
in  a  minimum  time  T  *  T. 
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Discussions  which  are  concerned  with  the  existence  of  the  opti¬ 
mal  controls  and  the  necessary  conditions  vfeich  the  optimal 
controls  must  satisfy  can  be  made  on  the  basis  of  the  problem 
specified  by  Equation  (1.9).  However,  to  gain  insight  from 
numerical  results  and  to  proceed  with  analytical  calculations, 
it  is  more  fruitful  to  work  with  the  simpler  models* 

This  study  has  three  main  objects.  The  first  is  to  formulate  a 
general  method  for  treating  the  many  possible  IMU  error  para¬ 
meter  configurations.  The  second  is  to  demonstrate,  through 
numerical  experimentation  with  realistic  error  parameter  con¬ 
figurations,  under  which  physical  conditions  it  is  possible  to 
identify  the  error  parameters  and  that  a  realistic  trajectory 
optimization  problem  does  exist.  The  third  object  is  to  formu¬ 
late  the  optimal  control  problem  and  discuss  the  existence  of 
the  optimal  controls,  the  necessary  conditions  which  the  optimal 
controls  must  satisfy,  and  the  computational  aspects  for  obtain¬ 
ing  the  optimal  controls. 
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CHAPTER  II 


GENERAL  APPROACH  TO  INFLIGHT  IDENTIFICATION 
OF  IMU  ERROR  PARAMETERS 


2.1  INTRODUCTION 

There  are  many  inertialiy  guided  vehicles  serving  as  a  common 
carrier  for  a  smaller  inertialiy  guided  vehicle.  Typically,  the 
smaller  vehicle  is  deployed  either  during  the  common  carrier’s 
travel  or  at  the  carrier's  nominal  terminus.  As  a  contrivance  to 
introduce  the  problem,  consider  the  Saturn  IV/ Command  and 
Seryice  Module  ( CSM  ).  Prior  to  launch  of  an  Apollo  mission  it 
may  be  difficult  to  calibrate  and/or  align  the  GSM  (slave)  IMU 
because  of  its  physical  arrangement.  The  master  IMU  ( IMU  used 
throughout  boost  powered  flight)  is  assumed  to  be  accurately 
aligned  and  calibrated.  When  the  vehicle  is  put  into  operation,  the 
slave  IMU  could  be  activated  prior  to  its  being  used  and  its  output 
observed.  The  observed  outputs  of  the  slave  IMU  and  the  master 
IMU  would  then  be  compared  and  the  difference  in  the  observations 
used  as  a  basis  for  inflight  calibration  and  alignment  of  the  slave 
IMU  with  respect  to  the  master  IMU. 


Under  ideal  conditions,  the  readings  from  the  two  IMU's  would  be 
the  same.  However,  there  are  several  factors  which  cause  the 
readings  between  the  master  and  slave  IMU's  to  be  different. 

These  factors  include  master  and  slave  IMU  error  parameters 
(gyro  drift  errors,  accelerometer  scale  factor  and  bias  errors, 
platform  and  accelerometer  misalignments),  accelerometer 
quantization  errors,  computer  errors,  and  random-induced 
accelerations.  Assuming  the  variances  of  the  master  IMU 
error  parameters  are  negligibly  small  so  that  the  master  IMU 
accelerometers  read  the  true  sensed  acceleration,  and  assuming 
availability  of  appropriate  probabilistic  descriptions  of  the  vehicle 
vibrations  and  of  the  computer  and  accelerometer  quantization 
errors,  which  are  random  in  nature,  various  techniques  may  be  used 
to  identify  the  IMU  error  parameters.  In  this  chapter  a  general 
approach  is  formulated  for  identifying  these  IMU  error  parameters  r 

2.2  Dynamic  Equations 

The  master  IMU,  denoted  by  IMU0 ,  is  the  reference  for  calibrating 
and  aligning  the  slaved  IMU,  XMUg .  It  is  assumed  that  EMUS  is,  in 
general,  a  strapped-down  system,  thereby  insuring  the  necessary 
generality  in  the  dynamical  equations  for  treating  either  gimbaied 
or  strapped-down  IMU's.  In  the  discussion  of  this  section,  the 
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notation  listed  below  is  used.  Later,  whenever  the  meaning  is 
clear,  a  less  complicated  notation  is  adopted. 

Referring  to  Figure  2.1,  we  let  denote  an  object  in  the  system 
with  orthogonal  axes  -1,2,3  ,  centered  at  B.P.^,  a 

point  fixed  in  the  object.  The  superscript  j  will  be  m  for  the 
vehicle,  o  for  1MU0 ,  s  for  £MUS ,  r  for  some  generally  non¬ 
stationary  reference  point,  and  i  for  an  inertially  fixed  reference 
point,  Jbp  is  the  vector  from  B.P.^  to  where  as 

above  j,  k  =  m,o,s,r,  and  i.  Letting  subscript  c  denote  a 
nominal,  or  assumed  position,  we  define  eb^i,s^  as  the  error 
vector  in  the  nominal  position  of  the  slaved  1MU,  and  ebf* » °)  as 
the  error  vector  in  the  nominal  position  of  the  master  IMU.  For 
compactness,  we  shall  use  a  dot  in  the  superscripts  to  indicate 
situations  where  j  can  either  be  s  or  o . 

The  position  and  velocity  of  IMU^.j ,  with  all  vectors  expressed  in 
vehicle  coordinates  Xq^,  are  given  respectively  by 

b<‘,-> .  ^i,-)  +  /■) 

t/* '  ^  =  6^1,  '1  + 


where  the  dot  over  a  vector  denotes  a  derivative  with  respect  to 


.-Ai!.- 


■Wl".  1  l.ll 


Figure  2. 2  Rotations  Between  Coordinate  Systems  Using 
Body  Angles 
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BMW 


inertial  space  and  the  small  circle  over  a  vector  denotes  a  deri¬ 
vative  with  respect  to  vehicle  coordinates,  X^.  Thus  the 
acceleration  is 

fif**  =  bf*’  +  n(t)('^  +  wxft(t)M  +  ^xr^  +  wxf^  +  wxwxrf  ^ 

so  that  the  acceleration  of  the  points  p  and  p’  (Figure2.1)  with 
respect  to  inertial  space  are 

SLf  =  l!^  ^ +  M(t)^ +  +(wxwx+ wx)eb^  ‘K  n(t)(') 

where  the  subscript  (•)  denotes  either  p  or  p\  The  acceleration 
of  IMUS  with  respect  IMU0 ,  &af°  *B\  is 

4a(°’  s)=fefo»  s)+ H(t/°»  s)+  2wxA(t/°»  s)+(o?xu>x  +  wx)n(t]l0’ s) 

\  (os)  (2.1) 
+  (wx(J)x+wx)eb^,  ' 


Equation  (2. 1)  gives  the  accelerations  between  the  master  and  slave 
EMU’s  due  to  random  vehicle  vibrations  n(t),  and  the  errors  in 
knowledge  of  the  nominal  locations  of  the  XMU’s.  We  can  rewrite 
(2.1)  as  Equation  (2.2)  using  the  following  notation: 


l«i  = 


ro 

U)g  0  *“  W  J 
■«2  U1  0 . 
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in’!  =  \[(o]  feu]  +  fcujj 

mi  =  [i  :2  [u\  *;  imj 

l(t)  =  (1  (t) ,  n  (t) ,  nj* 

aJ°,S)  *  mi  n(t)  +  in’]  e]j£’s)  +  i°,s) 

(2.2) 

Thus  the  difference  in  acceleration!?  between  the  master  and  slave 
IMU's  is  composed  of  a  term  due  to  random  vibrations  of  the 
vehicle,  tn]^(t) ;  a  term  due  to  an  error  in  knowing  the  true 
static  locations  of  the  IMIT’s,  [fl‘l  ejbf°,s^;  and  a  term  due  to 
gravity  and  the  angular  velocity  of  the  vehicle,  s^,  which  is 
known  as  a  function  of  time.  Next  the  term  is  considered 

in  the  context  of  the  measured  acceleration  errors. 

(i  •) 

The  acceleration  of  IMU^.j  with  respect  to  inertial  space  a'  *  , 

may  be  written  a 3  j|^  ^  =  j^  ;  +  a^  where  is  the 

(.) 

gravitational  acceleration  acting  on  IMU( . -j ,  and  a' 7  is  the 
thrust  acceleration  acting  at  IMU^.j .  Furihernore,  may 
be  written  as  =  4neasured  “  Aaf 80  we  can  write  the 
difference  in  acceleration  between  the  two  IMU’s  as 

-measured  ”  “  ( ^-measured  "  A-ineasured^ 

+  in]  n(t)  f  In']  +  j^0*8) 
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Henceforth,  we  shall  use  the  subscript  m  for  the  word  measured. 
We  note  that  the  vector  a{°>  8  ^  is  the  quantify  which  may  be 
physically  observed,  although  perhaps  indirectly,  through  the 
integral  of  af^  s^.  We  can  now  write  the  differential  equation  for 
the  measured  accelerations. 

The  motion  of  the  vehicle  is  described  by  the  equation 

E  -  g(E,%)  +  1(H,R,Ls) 

where  R  ,  R ,  R ,  represent  the  position,  velocity  and  acceleration 

of  016  vehicle  relatlTe  *°  “  inertial  reference  frame’  ■%)  48 
the  acceleration  due  to  gravity,  J(R ,  R ,  lg)  is  the  true  sensed 

acceleration,  Eg  is  a  vector  of  gravity  error  parameters,  and  t g 

is  a  vector  of  sensed  acceleration  error  parameters. 

M 

If  R'  '  is  the  true  position  of  the  geometrical  center  of  IMU^  then 
r  «  -  R^0)  satisfies  the  differential  equation 

d2r/ dt2  =  (£<s>  -  V  (is)  -  i0)), 

and,  neglecting  the  effect  of  the  error  parameters 

d2i/dt2  =  -  jf0))  +  a^»s^+  (Aa^  -  oa^) 

or 

d2r/dt2  =  fnj  n(t)  +  (Q’J  (2.3) 
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If  r(')'  is  the  computed  or  measured  position  of  IMU^j  and 
X*  -  -  B^>  then  Equation  (2.3)  may  be  rewritten  in  terms 


of  r*  as 


that  is 


d2ryat2=  tnil(t)  +  (n'J  s)  - s>  +  ^°> s)’  -  gf°» s) 

(2.4) 

If  Xq  =  b^°’  8\  the  nominal  or  reference  distance  from  IMU0  to 
IMUS ,  the  term  b^°*s)  can  be  written  as 

Jio  *  -  l(R(o0))]  +  wx^xii, 

It  is  reasonable  to  assume  that  g(B^  -B^)  s  j£(B^  -  B^)> 
so  that  Equation  (2.4)  becomes 


4n’  ^  =  (n}  2(0  +  fjyi  e^o  +  M  £q  +  A-m  S)  (2‘ 5) 

Further  we  can  define  the  observable  quantities  and 

B\  so  we  can  write  Equation 

(2.5)  a* 
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Next  it  is  shown  that  the  measurement  errors  can  be  written  as 
linear  combination  of  constants:  Aa^  =  B(t)^ ^ where  B(t)^  ^ 
is  a  matrix  relating  the  error  coefficients  of  IMUf .  j  to  trajectory 
variables  (acceleration,  angular  velocities,  and  time-varying  co¬ 
ordinate  transformations). 


Let  T^’  )  denote  the  transformation  from  an  inertial  coordinate 
system  to  a  moving  (or  body)  coordinate  system.  The  inverse 
transformation  is  written  as  T^*’  where  the  dot  as  before,  in¬ 
dicates  a  particular  object  in  the  system.  The  coordinate  systems 
specified  earlier,  xi*  \  stre  related  to  one  another  by  the  transformation 

2^'  =  where  a>0=  l>2>*  and  j,k  =  n,o,s,r,i. 

The  are  elements  of  the  usual  transformation  matrix,  which 

is  implicitly  defined  in  Figure  2.2  in  terms  of  vehicle  body  angles. 
Figure  2.3  is  included  to  provide  a  set  of  accelerometer  and  gyro 
coordinates.  In  addition,  two  transformations,  and  m(®),  are 
defined  to  express  the  transformations  of  the  master  IMU  and  the 
slave  IMU,  respectively,  relative  to  some  fixed  inertial  coordinate 
systems,  and  are  constant  transformations  relating  the  initial 
orientations  of  two  IMU's  to  some  inertially  fixed  coordinate  system. 


It  is  next  necessary  to  define  the  time-varying  transformation  M^(t) 
relating  the  slaved  IMU’s  present  position  relative  to  its  original 
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if  a  torque  is  applied  to  the 
inner  gimbal  about  the  x-axis, 
tht  wheel  and  inner  gimbal 
will  process  about  the  z-axis. 
If  a  torque  is  applied  to  the 
outer  gimbal  about  the  z-axis, 
the  wheel  and  inner  gimbal 
will  process  about  the  x-axis. 


Figure  2. 3  Accelerometer  and  Gyro  Coordinates 
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position  as  specified  by  (This  is  due  to  the  assumption  of  a 

strapped-down  system).  Similarly,  if  the  master  BSU  is  not 
inertially  fixed,  a  time -varying  transformation,  M^°\t),  is 
specified  in  terms  of  the  body  angles  ,  a  -  1 , 2 , 3 ,  the 

rotations  of  the  master  IMU’s  relative  orientations,  as  specified 
by  Then,  T^»')~  M^M(t)^,  and  a  vector  in  inertial 

coordinates,  v  j ,  can  be  expressed  as  a  vector  in  moving  (or  body) 
coordinates  v£  by  the  transformation  v^  =  vy).  Thus,  we 

may  write  the  thrust  acceleration  as  a£  -  AaB  =  (T  -  AT)(aj  -Aaj) 
where  the  superscripts  have  been  temporarily  dropped  for  ease  of 
notation,  a£  is  the  true  thrust  acceleration  in  body  coordinates, 

Aa£  is  the  error  in  measuring  the  thrust  acceleration,  expressed 
in  body  coordinates,  and  AT  is  the  error  in  knowing  the  true  trans¬ 
formation  T. 

Expanding  the  above  equation  and  then  neglecting  the  second-order 
terms  gives 

Aa£  -  a^  -  T^  +  TASj  +  ATaj 
or 

Aax  =  T_1AaB  -  T^ATaj  (2.7) 

Noting  that  the  term  ATT”1  can  be  written  as  an  antisymmetric 
matrix  of  misalignment  angles  expressed  in  body  coordinates, 
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ATT'1  =  [*B]  = 


0  -*z  *y 
*Z  0  -*K 

~*y  *x  3 


’  *B  =[*x  *y  *zj? 


we  have  from  Equation  (2. 7) 


-1 


A4i  =  T  (AaB  -  mjj)  ajj) 


Now  the  transformation  T  is  computed  by  integrating  the  equation 
d(T  -  AT)/dt  =  (T  -  AT)  (2. 8) 

where 


H- 


0 

"y 

<*>z 

0 

-Wx 

-Wy 

«x 

0 

jj 

and  wig  =  [wxjWy^zl 


is  the  nominal  angular  velocity  of  the  body,  and  where 


N- 


0  £y 

Ez  0  **ex 


-ey  ex  0 


,  eB  *  [Sjj,  ey,ezf  is  a  vector  of 


gyro  drift  rates  expressed  in  body  coordinates.  Expanding  (2. 8) 
and  neglecting  second-order  terms  gives 
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d(AT)/dt  =  CeB)  T  -  fufe]  AT 

-I 

but  from  the  definition  of  ATT  we  can  write 

(dT/dt)  f*j]  +  T  (dfyj/dt)  =  [Gg}  T  -  fo>B)  T  [*j] 


Next,  we  note  that  (dT/dt)  te|]  =  -fwg}.  T  }  which  implies 
T(df*jJ/dt)  =  feB)  T  or  dNtf/dt  =  T-1  [tB]  T  =  r^j 
Using  vector  notation,  this  latter  equation  becomes 
df^/dt  =  T  1  eR  =  £j 
whose  solution  is 


/t  i 

IW‘  eB(a)dc r  (2.9) 

*0 

This  means  Equation  (2. 7)  may  be  written  as 

-1  -1 
Aaj  -  T  (AaB  -  T  (#jJ  T  ag) 

or,  by  defining  an  antisymmetric  matrix  [aj]  whose  elements  are 
the  components  of  aj  ,  we  can  write 


-i 

Aaj  =  T  A£b  +  tejj 
Substituting  (2.9)  into  (2. 10)  gives 

AS!  =  T’1AaB+  (ajJ  ±  (to^  + 


(2. 10) 


T‘V)  £B(a)  da 


(2.11) 
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From  equations  (2.22)  and  (2.23),  infra,  we  have 


Ajin  =  B  (t)  K  and  -  B  Jt)  K 
a  —a  B  g  -=g 


where  Ka  is  a  constant  vector  of  accelerometer  errors  and  Kg  is  a 
constant  vector  of  gyro  errors. 


By  defining  a  vector  K  =  [^(t0)j ,  Ka,  Kg]*  and  a  matrix 

B(t)  =  Jtaj'j  T_1(t)  BgM  do  j 

we  can  write  (2. 11)  as  Aaj  =  B(t)K.  Thus,  for  the  two  IMU's 

Aa^  =  B^‘\t)|^^  (2.12) 

where  Aj|*^  is  a  vector  erf  errors  in  thrust  acceleration  due  to  the 
IMU{.)  instrument  errors. 

2.3  IMU  Instrument  Errors 

ffl 

The  accelerometer  measurement  error  for  the  i  accelerometer  may 
be  written  as 

iami  “  vi  +  [Ka’ °N- “Cj1  +  ~T  (Kgi'  %■  <2'!3> 

and  the  gyro  dr-t  rate  for  the  1th  gyro  may  be  written  as 
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(2. 14) 


£i  =  S  +  (K“l  \  V"*  lKI{  \  01  % 

+  *X2,j  a  *4  %  >  % 

In  these  equations,  a^  is  the  thrust  acceleration  in  the  ith 
accelerometer  (gyro)  coordinates;  is  the  bias  error  for  the 
i**1  accelerometer;  Kj^  is  the  scale  factor  in  the  ith  accelerometer; 

is  the  misalignment  of  the  i**1  accelerometer's  input  axis  in  the 
plane  of  the  input  and  normal  axes;  is  the  misalignment  of  the 
ifcil  accelerometer’s  input  axis  in  the  plane  of  the  input  and  cross 
axes;  w  is  the  angular  velocity  in  the  i  gyro  coordinates  (for  ideal 
gimbaled  systems  this  term  will  be  the  zero  vector);  is  the  1th 
gyro’s  constant  drift  rate;  Kw_  is  the  scale  factor  of  the  1th  gyro; 

YS.  (y  )  is  the  misalignment  of  the  ith  gyro’s  input  axis  in  the  plane 
defined  by  the  input  and  spin  (output)  axes;  Kj^  (Kg^ )  is  the  1th 
gyro’s  error  coefficient  due  to  mass  unbalance  along  the  input  (spin) 
axis;  ^Kg.  j  is  a  symmetric  matrix  where  the  elements  are  the 
1th  accelerometers  non-linearity  for  accelerations  along  the  itJl  axis 
and  the  elements  are  the  cross -axis  sensitivity  of  the  i**1  acceler¬ 


ometer  in  the  j  -  k  plane 


■  K  ^ 
- 


is  a  symmetric  matrix  where  K  *  j 


are  the  error  coefficients  sensitive  to  the  angular  velocity  com¬ 


ponents  w  •  for  the  i  gyro;  and  JKgJ  is  a  symmetric  matrix 
whose  components  Kg^  are  the  error  coefficients  sensitive  to  the 
acceleration  components  a  j  a^  (anisolastic  effect)  for  the  i^1  gyro. 
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It  is  certain  that  many  of  the  elements  in  the  above  matrices  will 
be  so  small  that  they  may  be  neglected.  If  the  gyros  are  used  on 
a  gimbaled  platform,  the  terms  proportional  to  u>  and  uP“  would 
be  neglected.  By  inserting  the  superscripts  o  and  s  in  (2. 13) 
and  (2.14)  the  corresponding  equations  are  obtained  for  the  master 
and  slave  IMU's,  respectively. 


In  the  above  error  equations,  the  input  accelerations  and  angular 
velocities  are  assumed  to  have  been  transformed  into  the  co¬ 
ordinates  of  the  particular  instrument  in  question.  Thus  assuming 
a  vector  vjj^  is  expressed  in  some  inertial  frame,  it  may  be 
expressed  in  accelerometer  coordinates  by 


x<*>  =  -  SW 

acc  It  i 

M 

or  in  gyro  coordinates  by  replacing  A'j '  by  G  j  ,  so  that 


) 


„(•)  .G(;)T 

gyro 

Applying  these  transformations  to  (2. 13)  and  (2. 14)  gives 


(2.15) 


and 


®i i->  ■  eoj)+!Kit>’  ysj»  4;’.  °i  4()*(t 


(2.18) 


*  4>*  4()*<')’  ^  ()-T 


3  8 


when  a^  is  the  thrust  acceleration  in  some  inertial  reference 
frame  and  ^00  is  the  an§:ular  velocity  in  the  same  inertial  reference 
frame. 


To  express  the  acceleration  measurement  errors  and  gyro  drift 
errors  in  platform  coordinates,  the  components  Aamj,  and  ej  are 
multiplied  by  the  transformation  4 1  and  G  /  respectively,  thus 


where  the  bar  under  A ,  G  indicates  that  only  the  appropriate  elements 
of  and  G^  are  included  here. 


If  (2. 15)  is  exuded,  it  could  be  written  as  the  scalar  product  of  a 
(•) 

vector,  pX  ,  which  depends  only  on  the  nominal  trajectory  parameters 
**i 

and  a  vector  of  the  i**1  accelerometer  error  parameters,  i.  e. , 


Aam  “ 

mj  a* 1 


~ai 


(2. 19) 


Similarly, 


W  -  41  4} 


(2.20) 
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*Pse*t}  WMN^ifSgS 


Thus  (2. 17)  and  {2. 18)  may  be  written,  respectively,  as 
aa(-)  =  [B<‘>(t)]K<‘> 


and 


~a 


(2.21) 


(*)  (*)  (•) 
gV  '  =  ! B g  (t)]  Kg 


Letting 


(2.22) 


where 


K(  )  _  fK(  )  M  K(  )1  K(  )  _  fK0  K( )  K(  )> 

-a  -  &g  ~  [ag^gg^gg] 

«s>  •  •  •  •]■ 

*8  ■  pv"1-!  -’8  •  •  ■  •K<i’  ■  •  ■]' 


(2,23) 


(2.18),  (2.19),  and (2.20). 


and  the  matrices  CB”(t)],  (Bg(t)]  are  defined  by  equations  (2.17), 


In  any  realistic  problem,  a  good  many  the  components  hi  (2. 13) 
and  (2, 14)  would  be  zero,  especially  for  the  master  IMU  and  for 
quadratic  error  terms.  Equations  (2. 21)  and  (2.22)  are  a  useful 
form  of  the  instrument  error  terms  for  use  in  the  minimum- 
variance  estimator  since  the  vectors  are  constant 

vectors  depending  only  on  the  error  parameters  c'  the  accelerometers 
and  gyros. 
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2. 4  Structure  of  the  Minimum-Variaice  Estimator 


It  is  necessary  to  assume*  that  a  statistical  description  of  random 
inputs  due  to  vehicle  induced,  random  acceleration  is  aval1  able,  and 
that  these  inputs  can  be  represented  by  white  noise  passing  through 
a  linear  dynamical  system.  Such  a  linear  filtering  operation  is 
performed  by  what  is  usually  called  a  shaping  filter  (Appendix  "A" 
presents  the  method  which  the  random  vibrations,  jj(t),  can  be 
represented  by  a  linear  dynamical  system  along  with  illustrative 
examples).  Thus,  it  is  assumed  that  the  vehicle  induced,  random 
accelerations  may  be  represented  by 

djj'/dt  =  Ap(t)n/  +  Bp{t)  w<t)  (2.25) 

where  Ap(t)  and  Bp(t)  are,  in  general,  time  varying  matrices 
which  are  related  to  the  covariance  matrix  of  n ,  and  w(t)  is  a 
white  noise  vector  process.  The  prime  on  n  is  to  indicate  that 
v’  may  be  composed  of  v  and  other  additional  suitably  defined 
vectors  arranged  one  upon  the  other.  We  shall  write  n  ■= 
to  indicate  that  only  the  components  n  are  used  from  v' .  We  now 
combine  equations  (2.6),  (2.12),  and  (2.25)  and  get 


The  observations,  or  measurements,  are  assumed  to  be 
and  (Actually  one  might  also  include  the  vehicle's 

angular  velocity  as  an  observation,  but  one  can  show  that  the 
angular  velocity  equation  can  be  uncoupled  from  the  rest;  thus, 
it  can  be  studied  independantly.)  These  observations  are  assumed 
to  have  errors,  and  it  is  assumed  that  the  error  in  the  observa¬ 
tions,  because  of  observation  error  parameters,  may  be  written 
as 

Az  =  Bm(t)km  (2.27) 

where  az  is  a  vector  denoting  the  error  in  the  observation  vector, 
z ,  due  to  the  observation  error  parameters,  km,  which  are 
assumed  to  be  constants  over  a  particular  observation  period,  and 
Bm(t)  is  a  matrix  relating  the  observation  error  parameters  to 
the  nominal  (or  known)  variables  in  the  observation  and  dynamical 
equations.  The  vector  km  would  include  such  things  as  biases  and 
scale  factor  errors  in  the  various  EMU  integrators. 

Next  the  observations  are  assigned  to  be  contaminated  by  a  random 
noise  v.  To  be  general,  it  is  assumed  that  this  noise  could  be 
correlated,  so  it  is  expressed  as 

v  =  +  vw  (2.28) 

where  is  a  vector  of  correlated  noise  components  and  vw  is  a 
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vector  of  white  noise  components  ^  As  above  (2.25),  it  is  assumed 
the  correlated  noise  may  be  represented  in  the  Markovian  form 


dj'J./dt  =  Am(t)j£  +  Bm(t)v^,  (2.29) 

where  the  prime  carries  the  same  connotation  as  above.  The 
observations  are  thus 

+  Bm(t)km  +  (rlj£  +  Vw  (2.30) 

Defining  the  state  vector  as 

Equations  (2.26)  and  (2.30)  may  be  written  in  canonical  form  as 
X  =  AX  +  BW  +  1 

Z  =  HX  +  V  (2.31) 


where 
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Equations (2.31)  have  the  discrete  form  X(k+1)  =  <?(k+l,k)X(k)  +  W(k) 
and  Z(k+ 1 )  =  H(k+ 1  )X(k+ 1 )  +  V(k+ 1 ) .  We  next  show  how  the  dis¬ 
crete  equation  set  is  formulated  in  terms  of  a  specific  problem. 

.5  Specialized  Calibration  and  Alignment  Problem 

In  this  concluding  section  a  discrete  form  of  (2.31)  is  developed 
in  a  manner  to  give  concrete  meaning  to  the  ideas  presented  above, 
and  to  indicate  the  approach  used  in  the  next  chapters  where  further 
examples  and  numerical  results  are  presented.  In  the  problem 
considered  here,  IMIT0  is  well  calibrated  and  aligned  prior  to 
vehicle  operation  so  that  IMUS  could  be  considered  adequately 
calibrated  and  aligned  when  compared  with  IMU0  as  a  standard. 
Thus,  the  vector  of  IMU0  error  parameters,  k(°),  may  be  taken 
to  be  the  zero  vector. 


Next  it  is  assumed  nonlinear  error  terras  and  fixed  misalignment 
angles  on  the  accelerometers  and  gyros  are  negligible  and  that 
these  instruments  are  nominally  located  so  that  the  instruments’ 
sensitive  axes  form  orthogonal  triads...  Under  these  assumptions 
the  matrix  3^s^(t)  becomes 


B(S)(0 


(a]!M^S)(t)  M(S)^ 


da  )  M^s)$ 
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a  3x21  matrix  where 


A  = 


1 

0 

0 


0  0  ax  0  0 

1  0  0  ay  0 

0  1  0  0  az 


9  = 


l 

o 

0 


0 

1 

0 


0  Ujj  0  0  ax  2y  0  0  0  0 

0  0  Uy  0  0  0  iy  o  0 

1  0  0  W2  0  0  0  0  ag  ax 


and  kfs^ 


is  given  by  the  21 -vector 


,  iS),  J&>>  £(oS>, 


In  the  matrix  B^s)(t),  ax,  Cy,  a2,  u^,  wy,  wz  are  nominal 
acceleration  and  angular  velocity  components  in  the  IMUS  co¬ 
ordinates.  Assuming  the  random  angular  acceleration  of  the 
vehicle  is  negligible,  u>  small  or  nearly  constant,  means  the 
"angular  velocity  matrix"  is  defined  by  [n'3=  [till  2  [wj],  and  by 
drfining  a  vehicle  vibration  vector  as  v(tf  =  then 

(2.6)  becomes 

j.  l£ml 1 1  MX.  -°. _ i 

dt  UfraJ  L-*ain+2+Hi’ 


Next  it  is  assumed  n  is  suitably  approximated  by 

dn/dt  =  [Ail.  Jl  >12]  n  +  B*(t)  w 
L  A21  ,  a22J  ~ 


where  n*  =(.22,213]  where  2.2  and  £3  are  3-vectors  and  the  Ajj  are 
3x3  time  varying  matrices,  ng  =  I,  Jlj  =  fi  -  B’(t)  is  a 

6x6  time  varying  matrix  and  w  is  a  6 -white  noise  vector.  Thus 
we  can  write  j?  =  Ap(t)  +  Bp(t)  sl  where  ZL'*=  (nj,  ng, 
and 


AP(t)  * 


[rJ  1  0 


An : a,2  ,  Ep(t)  *  ...  irwinlo] 

*  t  *  *3  (tj 

A21  *  A22J  L  6  I 


where  fl]  is  a  6x6  identity  matrix,  m  is  uncorrelated,  and  kfm^=0. 
Assuming  the  angular  velocity  equation  is  treated  separately,  the 
observation  equation  becomes  Z  •-  t]  +  v(t)  where  v(t)  is  a 
6-vector  of  white  measurement  noise.  If  we  define  a  36-state 
vector  as  X*  =  v^,  ?£,  fL>  H3,  and  the  matrices 


A(t)  =  i- 


0 

I 

0 

0 

0 

0 

rr] 

-b(s) 

0 

0 

AP 

0 

0 

0 

0 

0 

,  B(t)»  BP 

:o: 


and  H  =  |(T]jOj,then  the  continuous  problem  is  specified  by 
X  =  A(t)  X  *  B(t)  w  and  Z  =  HX  +  v(t).  The  equations  for  the 
discrete  formulation  are 


■  M.' 


Xk+1  =o(k+l,k.)Xfc  +  Sfe  and  £k+i  =  H(k+2)X(k+l)+yk+i 

(2.32) 
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where  <t»(k+  l,k)  is  given  by 


■^k+i,k) 

♦Jk-fl.k) 

^k+l,k)‘ 

<t>(k+  l,k)  = 

0 

♦4(k+l,k) 

0 

0 

0 

II  . 

(2.33) 


The  submatrices  of  (2.33)  are  given  by 


<^k+l,k)  = 


I 

!  AM 

"  "0" 

;  ’1  \ 

<4(k+l,k)  = 


m  m 

Atl 

I 

-k+1 


^k+l,k)  =  f+[(<T,k)M{mk<T%U 
/-k+1 

4>3(k+l,k)  ^♦•(a.k)B(8)(cr)d(7 

1 1  is  a  21x21  identity  matrix  and  ^(k+  l,k)  is  the  solution  of  the 
homogeneous  equation  d  ^(t,  r)  / dt  =  A^(t)  ^(t,r)  where,  as  in¬ 
dicated  below,  ^(t,t)  =  I.  Qualitatively,  ^  is  a  matrix  which 
gives  the  contributions  to  and  vm  due  to  vehicle  vibrations  l , 
and  \  is  the  9x9  transition  matrix  for  the  random  vehicle 
vibrations.  The  white  noise  sequence  w^  is  given  by 

=  /k+1  £  «4((t,  k)B^r)  x(oj|  dcr  .  The  wliite  noise  observation 
sequence  y^j,  is  similarly  obtained  from  the  whilw  noise 
observation  process,  v(t).  Ij  =  identity  indicates  that  the  error 


parameters  remain  constant.  The  transition  matrices  have 
the  following  properties: 

1)  ^(t3,  t2)^(t2,  tx)  =^(t3,  tx)  for  all  tltt2,t3, 

3)  ^(t ,  r)  =  X(t)X(r)_1  where  X(t)  is  the  solution 
to  the  equation  dX/ tit  =  A(t)X(t),  X(o)  =  I. 


From  the  above  properties  it  is  seen  that  if  to=0,  *(t1,0)  =  X(o), 
and  it  is  not  difficult  to  see  that  in  general  t^)  =  x(k+l)  (tj) 

"  ^  -  <tj,  0)  where  X^  \t^)  is  the  solution  to  the  matrix 

equation  X(t)  «A(Ut^)X(t)y  X(o)  =1.  This  means  that  it  is 
only  necessary  to  solve  for  the  matrix  X'tj)  for  each  time 
interval,  (tk,  t^),  with  an  updated  A(t)  matrix  using  the  same 
initial  conditions  X(o)  =  I.  Thus  the  computer  routine  for 
generating  the  transition  matrix,  <Ktfc+2>  t^)  is  not  changed  for 
each  k,  only  the  elements  of  A{t)  need  to  be  updated. 
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CHAPTER  HI 


FORMULATION  OF  SPECIFIC  ERROR  PARAMETER 
IDENTIFICATION  PROBLEMS 

In  Chapter  I ,  a  very  general  formulation  of  the  IMU  error  parameters 
was  considered.  This  formulation  is  general  in  the  sense  that  most 
specific  problems  can  be  placed  into  this  framework.  When  a  specific 
problem  is  considered,  the  equations  of  Chapter  n  simplify  considerably 
and  the  dimensions  of  the  various  matrices  are  decreased  accordingly. 
For  example,  it  could  be  assumed  that  all  the  error  parameters  are 
known  and  that  only  the  initial  misalignment  angles,  $(t0),  need  to  be 
estimated.  Another  important  simple  example  might  be  to  e«*timai.i 
the  random  motion  only.  This  might  occur  when  a  vehicle  is  launched 
from  the  wing  of  an  airplane  or  mother -ship. 

In  this  chapter  we  discuss  three  important  IMU  error  parameter  con¬ 
figurations  which  were  considered  in  this  study.  They  are  as  follows: 

1.  Initial  misalignment  angles,  ± ,  only. 

2.  Initial  misalignment  angles  plus  gyro  mass -unbalance 
terms,  ks  and  kj  . 

3.  Initial  misalignment  angles  plus  gyro  constant  drift 
rates,  e . 
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Both  time  correlated  and  uncorrelated  random  accelerations  are  con¬ 
sidered.  The  correlated  model  used  is  discussed  in  considerable  detail 
in  Appendix  "A”.  The  correlated  random  acceleration  is  assumed  to 
act  along  one  axis  of  the  vehicle  (axis  of  major  disturbance)  and  depends 
upon  specific  parameters  which  cannot  be  readily  varied  for  parametric 
or  analytic  studies.  For  this  reason  extensive  experimental  results  for 
this  model  are  not  included  -  although  they  are  available.  The  main 
purpose  of  this  study  is  to  determine  effects  of  state  and  observation 
noise  levels,  and  the  trajectory  parameters  on  the  identification  of  the 
errors,  and  the  uncorrelated  noise  models  are  most  convenient  for  these 
purposes. 

For  these  specific  models,  the  master  IMU  is  assumed  to  be  perfect, 
and  the  slaVved  TMU  is  assumed  to  be  a  gimballed  system.  The  only 
degree  of  freedom  in  the  nominal  trajectory  specification  is  assumed 
to  be  in  the  vehicle's  pitch  angle.  That  is,  the  vehicle  is  assumed  to 
have  a  fixed  thrust  program.  The  experimental  results  of  this 
chapter  indicate  that  the  quality  of  the  identification  depends  very  much 
on  the  manner  in  which  the  particular  error  parameters  enter  the  sys¬ 
tem  and  on  the  power  levels  in  the  random  disturbances.  Thus  a  mean¬ 
ingful  optimization  problem  can  be  formulated  with  respect  to  the  nominal 
trajectory  variables,  such  as  for  example,  the  trajectory  pitch  angle. 
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3. 1  Time  Correlated  Random  Vehicle  Accelerations 

For  the  present  fur  poses  it  is  assumed  that  significant  random 
motion  occurs  along  the  vehicle's  yaw  axis  only.  Such  a  random 
motion  might  occur  if  the  thrust  engine  were  to  move  randomly 
about  the  pitch  axis  while  controlling  che  vehicle's  trajectory  to  lie 
in  a  vertical  plane.  Furthermore,  it  is  assumed  that  this  random 
motion  is  stationary  and  has  a  power  spectrum  of  the  form  discussed 
in  Appendix  "A",  Section  A.  l  .  In  this  section  the  assumption  that 
the  given  power  spectrum  may  be  suitable  approximated  by  the  fourth- 
order  rational  power  spectrum 


SM 


V>  +  *pb 

(a>  -  a)^  +  b^  (w  h  a)^  + 


(3.1a) 


is  made.  The  corresponding  autocorrelation  function  is 

Rx(r)  =  *0  exp(-bM)  cos  aT  (3.1b) 

By  Equation  (A.  5)  of  Appendix  A,  a  model  of  the  process  defined  by 
Equation  (3.1a)  is  given  by 


The  "noise",  w(t),  is  such  that  E[w(t)]  =  E[w(t)  w(r)j  =  j(t  -  t).  j 

For  engineering  purposes,  Equation  (3.1c)  is  usually  cliscreti2ed  by  j 
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approximating  the  white  noise  process  w(t)  by  the  white  noise 
sequence  / yftT  where  k  denotes  the  number  of  the  sample  and 
a  denotes  the  sampling  interval  {see  Sec'lon  A. 4).  The  discrete 
form  of  (3.1c)  is  given  by  (A.  20c)  as 


Simplifying  the  notation,  Equation  (3. Id)  may  be  written  as 

Xi  =  An  *i  -  1  +  bn  wi  _  i  ( j.  1) 


where  the  definition  of  Xj,  Ag,  and  b jj  is  obvious  from  Equation 
(3.1c).  The  ’’white  noise”  sequence  wj  _  j  has  the  property 

E  [wj]  =  0  and  cov{Wj,Wj)=  E  [(Wj  wj)]  =  (3.1e) 

Another  discrete  model  of  the  process  defined  by  Equation  (3.1a) 
is  given  in  Section  A.  4,  Equation  (A.21c),  as 


To  simplify  the  notation,  Equation  (3.2a)  is  written 


(3.2a) 


*i  =  Ajjj  Xj  „  i  +  bnr  wj  _  i  (3.2b) 

The  definitions  of  ai*  a£,  bj,  bg  are  given  in  Section  A.  4. 
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Next  it  is  assumed  that  the  random  motion  which  is  introduced 
through  the  "angular  velocity  matrix"  [&*]  is  negligible.  This  is  a 
reasonable  assumption  if  the  nominal  angular  velocity  u  is  small 
and  nearly  constant.  In  the  present  case*  =  (O,^0’  0)* 


where  e 


(otm) 


is  the  vehicle  pitch  angle  and  §2 


(o,m) 


is  zero 


except  at  the  transitions  where  it  is  a  constant.  This  means  that 
the  "angular  velocity  matrix"  [n]  is  defined  by 


M  =  I 


Xl2  [w]I  * 
1  » 


1  0  0  !  0  -2wj  2u« 

1  * 

0  1  0  J  2o>3  0  -2^ 


1 ' 

1 

2yl 

0  j 

0 ! 

0 

0 

0 1 

! 

0 

0 

0 

l! 

"2w2 

0 

0 

^3. 3a) 


and  the  vehicle’s  random  acceleration  vector  is 

a(t)  >  (0  0  s  0  0  a]* 


(3.3b) 


According  to  Equation  (3.1c),  n  and  H  may  be  expressed  in  the  form 

/ 


_d_ 

dt 


n 

ft 


n 


all  a12  0 

a21  *22  0 
1  0  0 


\f 


f ? 


'll 


ft  J  +  j  b21  |  w  (3.3c) 

/  1  0 


where  n’  is  an  additional  variable  which  is  required  to  taJke  care  of 
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tile  correlation  in  1?.  The  elements  in  the  matrix  equation  above  are 
obtained  from  Equations  (3. 1c)  or  (3. 2a) ,  Equation  (3.3c)  may  be 
written  in  discrete  form  as 

ail  aI2  0  /  &  /  P  ll ' 

0/21  ®22  0  I  *  +  *21  wk  (3.3d) 

A  0  1  i  \  /k  0 

where  the  subscripts  k  +  1,  k  refer  to  the  sampling  interval,  and  the 
elements  in  the  matrices  are  obtained  from  Equation  (3.  Id)  or  (3 .2b). 


The  random  motion  of  the  vehicle,  in  vehicle  coordinates,  may  thus 
be  written  as 

n*m)  =  in)  \ti£  =[i|2Mj  [I  0}  (0  0  n  0  0  n  IPJ* 

If  we  substitute  into  the  above  matrices,  this  becomes 


(3.3e) 


(3.4a) 
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where 


*  0  0  202°,m) 

fO"}  s  0  0  0  .  and  «*  -  (K,  n«,  n  )*  (3.4b) 

.10  0 

• 

In  (31  it  is  shown  that  under  reasonable  assumptions  the  effects  of 
/  _  ^  \ 

@2  may  be  neglected;  that  is,  random  coriolis  accelerations 
can  be  neglected. 

3.2,  Uncorrelated  Random  Accelerations 

For  various  reasons,  it  may  be  possible  to  assume  that  the  discrete 
form  of  the  random  acceleration  due  to  vehicle  vibrations  may  be 
uncorrelated  in  time.  Such  would  be  the  case,  for  example,  if  the 
estimation  interval  A  was  required  to  be  greater  than,  say  0. 1  sec. 
Such  a  requirement  would  depend  on  how  fast  a  digital  computer 
could  process  the  data,  i.e. ,  the  velocity  differences  between  the  two 
IMU’s.  In  the  present  study,  it  has  been  found  that,  in  order  to 
estimate  initial  misalignments  and  mass -unbalance  drift  rates,  A 
would  need  to  be  greater  than  0. 1  sec.  In  the  event  that  A  >  0. 1 
Figure  3. 1  illustrates  how  the  correlated  noise  model  assumed 
in  the  previous  section  might  generate  essentially  uncorrelated 
velocity  sequences.  As  a  further  example,  a  sensible  candidate 


Figure  3. 1 


V 


COBKIA7ICM7MS,* 


Illustrations  of  the  Band-Pass  Correlated 
Noise  Model 
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Cci*  a  suboptimal  estimator  would  be  one  in  which  the  random 
accelerations  arc  assumed  uncorrelated  in  time.  It  is  thus  assumed 
that  the  random  acceleration  added  to  the  IMU  velocity  differences 
is  of  the  form 


rll  0 

0  r22 
0  0 


(3.5a) 


2 

Further,  it  is  assumed  that  aR  -  -  **22  ~  **33*  TW® 

assumption  is  made  for  computational  convenience  and  is  not  con¬ 
ceptual  in  nature.  The  random  acceleration  in  IMUg  coordinates 
is  given  by  |with  T  =  T^0>  T^m’  , 


£k 


(Sv  _  Ak+  1)A  (0>S)  (m>0) 


'■l 


r(m) 

(t)  dt  k 

Va- 


(3. 5b) 


The  covariance  of  this  random  acceleration  is  given  by. 


Cov(rjc,  rj)  =  E 


\  i 

ffi+i  \ 

=  E 

T(t)  r(t)  dtl.j 

I  T(t)  r(r)  drl 

ytk  1 

/ 

(3.  5c) 


due  to  the  fact  that  T(t)  is  an  orthogonal  transformation. 


3.3  Initial  Misalignment  Angles,  i  =  1,2,3 


If  it  is  assumed  that  the  accelerometer  error  parameters  and 
the  gyre  error  parameters  k£8^  are  either  zero  or  known  (Chapter 

o 

II),  then  the  difference  in  acceleration  between  the  slave  and  master 
IMU  (expressed  the  IMUS  coordinates)  becomes 


J8>0)  =  kmis  [a^  *  +  t 


In  this  equation  the  terms  are  defined  in  Chapter  II. 


(3.6a) 


If  the  random  coriolis  acceleration  is  neglected  (as  described  above) 


(3.6c) 


[n"l2'=io][S  (3.6b) 

In  order  to  write  Equation  (3.5b)  in  discrete  form,  let 

=  vf8’?^  (  Aafs,0^(7)  dr  (3.6c) 

~  Jk-1 

where  k  denotes  the  time  t  =  t^.  Equation  (3.6c)  may  then  be 
written  as 

4°’8)  =  +  kmls  t  [a(SV)j  dr*  ^T(m’  S)(r)drimVl  (3.6d) 

Combining  Equations  (3.6d)  and  (3.2b)  into  one  equation  gives 


e  i 
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an  observation  vector 


% 

s 

>1 


p 

Cf 

i 


yk  .  [y(°’s),  y<°-8)], 


(3.8b) 


a  state  transition  matrix 


4> 


k.k-l  & 


XV 

/  T(m* 8)  (t)  dr  [a»  ] 


k-l 


in 


kmi.  /  |»W(T)]dT 

k-l 


a  measurement  matrix 


M  -  tl  i  0!  ;  I  is  3x3  and  D  is  3x5 


(3.8c) 


(3.8d) 


and  a  process  noise  vector 


Ik  " 


’m 


Wt, 


(3.8e) 


With  these  definitions,  Equations  (3.6c)  and  (3.7a)  become 


-k  ~  *k,  k-l  -k-l  + 

Ik  “  M^k  + 


(3.9) 

(3. 10) 


with  the  covariance  matrices  of  wjj  and  r^  given  by  and 
respectively. 
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In  the  discussion  which  follows,  it  is  convenient  to  define  the 
quantities 

(k+  1)4/ 

AVi.m.(k+  ^  -  kmis^  ( V 
and 


)\dT* 


(3. 11a) 


(3.11b) 


In  this  section  we  obtain  the  equations  for  estimating  initial  gyro 
misalignment  angles,  ,  and  the  gyro  drift  terms  which  are  due 
to  gyro  mass  unbalances.  The  drift  rate  for  the  1th  gyro  is 


where  the  following  notation  is  used: 


(3. 12a) 


The  superscript  gj  refers  to  the  i**1  gyro; 

kj^  =  the  i^1  gyro's  error  coefficient  due  to  mass- 

unbalance  along  the  gyro's  input  axis,  in  deg/hr/g 

kg^  =  the  i**1  gyro's  error  coefficient  due  to  mass- 

unbalance  along  the  spin  axis,  in  deg/hr/g 

d&i)  -  the  sensed  acceleration  vector  in  ft/ sec5*,  expressed 

ih 

in  the  coordinates  of  the  i  gyro 
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^m.u.  =  0*15068493  x  10  ^  is  the  factor  which 
converts  deg/hr/g  to  rad/sec/(ft/sec^) 
e^i)  =  cirift  rate  of  the  gyro  in  rad/sec 

The  arrangement  of  the  gyro  coordinates  on  the  IMUS  is  shown 
in  Figure  3. 2 . 


After  introducing  the  required  coordinate  transformations,  the 
acceleration  error,  in  slaved  IMU  coordinates,  due  to  mass  un-- 
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The  total  velocity  error  due  to  mass  unbalance  drifts  between 
samples  is  then  given  by: 


The  matrix  product  in  Equation  (3.12c)  is 


'  °  'V«2 

a  a 
S2S3 

0  -a  a  a  a 

S3  s3  s2  S1 

r41j  •  [a^2  = 

W  .  ^  J 

° 

-a  a 
si  s3 

&  z  0  -a  a0 

s3  s2  S1  S1 

as2asx  asas2 

0 

~a  a  as  a 

S2  S2  ®1S3  0  J 

(3. 12d) 


where  the  superscript  s  has  been  dropped.  Thus  each  element  in 
the  integrated  3x6  matrix  of  Equation  (3. 12c)  will  involve  a  term 
of  the  form 


Ii>3(k+l)«-^+1^4S)(9|/taj(t')c!t’jdt;  i,j  =  1,2,3  (3.12e) 


It  is  convenient  to  define 


(k+  t)A 


. »]•/ 

kA 


|/  [aa(S)<T)]  ^  ^ 


(3  12f 
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so  that  Equation  (3.12c)  may  be  written  as 


a-v-m.u.<k+1<k>  ’  fcm 


.J4SW, 


m.u.—m.u. 


(3. 12g) 


The  estimator  equations  are  then  obtained  by  replacing 

““mis  s  J  mis 
in  Equation  (3.6c)  by 

[kmis  (v(sS>)mis  |  km.u.  m.u.] 

In  the  case  of  this  configuration ,  it  is  assumed  that  the  random 
accelerations  are  uncorrelated,  so  that  Am  =  0  in  Equation  (3.6e). 
Further,  the  vector  bjjj  in  Equation  (3.6c)  is  replaced  by  the 
3x3  matrix  of  Equation  (2. 5b). 


3.5  Initial  Misalignment  Angles  and  Constant  Gyro  Drifts 

In  this  section  we  obtain  the  equations  for  estimating  the  initial 
gyro  misalignments,  ,  and  the  constant  gyro  drift  terms 
,  i  =  1 , 2 , 3 .  The  constant  drift  of  the  i^1  gyro,  ,  is 
expressed  in  deg/hr.  The  arrangement  of  the  gyros  is  as 
shown  in  Figure  3.2.  The  velocity  error  expressed  in  IMUS 
coordinates  between  samples  will  thus  be 
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Avdr^l,k)  =  k d#r<jT  )A(iS)(t))  .  jjf  e  dt'J  dt 


=  k 


4Ti - 


(3. 13a) 


In  this  equation  M  is  the  matrix  defined  in  Chapter  II, 

-5 

Section  2.2,  and  =  0.48481368x10  is  the  conversion 

factor  from  deg/hr  to  rad/sec.  Next  we  define 


[4  ..£**■  fH 


dt 


(3. 12b) 


so  that 


AWk+1'k>  3  kd.r.  [Ts)]d.r.£ 


(3.12c) 


The  estimator  equations  are  then  obtained  by  replacing  ^mjS\vs  /m<g 
In  Equation  (3.6c)  by  [k^  I  r  (v^  J  .  For  this 

estimator  it  is  assumed  that  the  random  accelerations  are  uncorrelated, 
so  that  in  Equation  (3. 6e)  Am  »  0,  and  bjjj  is  replaced  by 
Equation  (3.5b). 


For  convenience  in  the  discussion  of  the  experimental  results,  the 
following  notation  is  used  for  the  configurations  of  the  various  slaved 
IMU  error  models : 
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Configurali.cn  XA:  Identification  of  initial  misalignment  angles 

±,  using  correlated  random  acceleration. 

Configuration  IB:  Identification  of  the  initial  misalignment 

angles  ±,  using  uncorrelated  random  accel¬ 
eration. 


Configuration  E:  Identification  of  the  initial  misalignment 

angles  * ,  and  the  constant  gyro  drifts  6  . 
Configuration  El:  Identification  of  the  initial  misalignment 

angles  ±,  and  the  mass -unbalance  gyro 


drifts,  k 

*  -m.u. 


i 
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CHAPTER  IV 

EXPERIMENTAL  RESULTS 

It  is  quite  useful  to  obtain  an  idea  of  the  type  of  results  which  are 
possible  for  the  trajectory  optimization  which  is  considered  in 
later  chapters.  In  this  chapter  we  give  some  experimental  results 
concerning  the  identification  of  the  error  parameters  for  config¬ 
urations  XA,  IB,  n  and  IH.  The  effects  of  state  noise,  measurement 
noise,  and  the  nominal  trajectory  are  considered.  The  nominal 
trajectory  is  reasonably  assumed  to  be  specified  by  the  thrust 
acceleration  and  by  the  trajectory  pitch  angle  (which  is  assumed 
very  nearly  equal  the  vehicle’s  pitch  angle).  The  vehicle  considered 
is  assumed  to  have  a  fixed  acceleration  profile,  so  that  the  only  degr  ee 
of  freedom  in  specifying  the  nominal  trajectory  is  its  pitch  angle,  n . 
These  functions,  m( t),  te  [tQ,  tQ+  T] ,  are  also  constrained  in  a 
certain  sense.  Generally,  the  pitch  angle  is  constrained  to  be  ±  10 
degrees  from  some  reference  value,  and  the  pitch  rate  might  be 
similarly  constrained.  Typically,  the  nominal  trajectories  would  be 
as  shown  in  Figure  4.1.  It  should  be  noted  that  most  o l  the  simplify¬ 
ing  assumptions  are  made  in  order  to  facilitate  the  parametric  study, 
and  arc-  not  conceptual  in  nature. 
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Figure  4-i  Reference  Thrust  Acceleration  Profile 
and  Trajectory  Pitch  Angle  Profile 


A  discussion  of  the  various  initial  conditions  (a  priori  estimates, 
etc.)  may  be  found  in  (3j.  The  notation  used  below  is  essen¬ 
tially  as  follows: 


4'1 

% 

ah 

W 


the  correlated  random  acceleration  power 

the  uncorrelated  random  acceleration  cov;triance  matrix. 

the  observation  noise  covariance  matrix, 

the  standard  deviation  in  the  estimated  misalignment  angles 

the  standard  deviation  in  the  estimated  constant  gyro  drift, 
rates. 

the  standard  deviation  in  the  estimated  mass -unbalance 
gyro  drift  terms. 
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The  initial  IMUS  orientation  angles,  do  not  seem  to  affect 

.  However,  the  of  the  individual  estimated  ^ ,  do 
depend  on  this  orientation.  For  this  reason  an  orientation  in  which 
the  accelerometers  each  sense  about  the  same  thrust  acceleration 
was  chosen  for  most  of  the  experimental  work  to  minimize  the  maxi¬ 
mum  of  the  o-^  .  The  geometry  is  as  shown  below . 


The  thrust  vector  in  IMU0  coordinates  is 


and  it  is  necessary  to  pick  the  three  fixed  rotations  4°,s^,  0^°’ 
0i}O,S^  suck  the  t^iree  direction  cosines  at,  otg,  agare  equal. 
Thus  , i  \ 


t(33,0’*>,  0J*0’*')  aj0* 
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If  4°’ m)  =  (TP  -90°)  =  -25°,  and  0^°’s)  =  TP  -45°,  the  solution  to  the 
above  equation  gives  0(°>s)=48. 9°,  e(°>s)=20°,  and  0|o>s)  =-22.5°. 


4.1  Experimental  Results  for  Correlated  Random  Acceleration 
(Configuration  IA) 

The  correlated  random  acceleration  model  is  studied  in  considerable 
detail  in  [3]  For  this  configuration  the  gyro  drift  rates  are 
assumed  negligible  so  that  only  the  initial  misalignment  angles 
are  estimated.  The  effects  of  the  sampling  rate  (data  processing 
rate) ,  A ,  are  shown  in  Figure  4.2.  The  sampling  time  is  varied 
between  0.025  sec  and  0. 100  sec  to  determine  how  much  improve¬ 
ment  in  the  identification  of  the  ^  and  the  correlated  noise  com¬ 
ponents  is  possible.  The  most  suitable  sampling  rate  depends  on 
the  particular  correlated  noise  model,  in  that  it  should  be  fast 
enough  to  obtain  a  reasonable  representation  of  the  correlated 
random  acceleration.  The  effects  of  the  random  acceleration 
rms  power  levels  and  measurement  noise  levels  are  shown  in 
Figure  4. 3  .  It  seems  that  it  is  of  advantage  to  structure  the 
estimator  for  correlated  random  accelerations  if  it  is  only  required 
to  estimate  the  initial  misalignment  angles  and  if  the  correlated 
noise  is  of  a  sufficiently  low  "frequency"  (say  less  than  three 
cycles/sec) .  It  is  also  of  advantage  to  consider  the  correlated 
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model  if  it  is  of  prime  importance  to  estimate  the  random  motion 
vdth  all  error  parameters  known,  as,  for  example,  in  launching 
from  a  mother  ship.  If  errors  in  addition  to  the  4^  are  also  to  be 
estimated  it  does  not  seem  possible  to  use  the  time  correlation  in 
the  random  acceleration  to  advantage. 

Since  the  results  for  the  correlated  noise  model  depend  on  many 
parameters  and  on  the  specific  noise  model  considered,  further 
discussion  will  not  be  made  here.  Considerable  experimental 
results  for  the  correlated  random  acceleration  model  may 


SAMPLING  TIME  A(SEC) 

Figure 4.2  Steady-State  Standard  Deviation  of  the  EsHmate 
of  the  Acceleration  Noise,  or  a  for  Various 
Sampling  Time  8  A 
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found  in  (  3  ].  In  this  case  a  parametric  study  is  not  as  readily 
made,  or  as  meaningful,  as  for  the  uncorrelated  noise  model 
which  is  considered  next.  Since  one  of  the  main  objects  of  this 
study  is  to  consider  the  maneuver  which  allows  that  the 
identification  of  the  errors  be  accomplished  in  an  optimal 
fashion  (Chapter  VI),  the  uncorrelated  model  is  used  because 
it  provides  a  better  setting  for  this  purpose. 


Figure  4.3a  Steady-State  Standard  Deviations  for  the 
Estimated  Random  Acceleration  n 
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STANDARD  DEVIATIONS  IN  THE  ESTIMATED  MISALIGNMENTS,  o-a  (ARC  MIN) 


4.2  Experimental  Results  for  Initial  Misalignment  Angles 
and  Uncorrelated  State  Noise  (Configuration  IB) 

In  this  case  the  data  the  data  processing  rate  was  fixed  at  A  =  0. 10 
seconds  and  the  ’’equal  acceleration”  orientation  discussed  above 
was  used.  The  effect  of  the  state  noise,  measurement  noise,  and 
trajectory  pitch  profile  are  considered.  Figure  4.5  illustrates 
the  sigma  4>;  for  different  values  of  measurement  noise  standard 
deviation,  The  random  acceleration  variance  o^,  for  these 

O 

cases  was  fixed  at  =  10"  ,  and  Nominal  Trajectory  IB-7,  Figure 

4.4,  was  used.  Figure  4.6  illustrates  the  sigma  ^  for  different 

2 

values  of  random  acceleration  variance,  cr^.  A  pertinent  discussion 
of  the  effects  of  process  noise  (random  acceleration)  and  measure¬ 
ment  noise  on  the  covariance  equation  P,  is  given  in  Section  4. 5. 

The  manner  in  which  the  IMUS  is  aligned  relative  to  IMU0  (through 
the  angles  8^°* and  the  trajectory  pitch  profile  enter  the  variance 
equation  are  outlined  in  Chapter  V.  Since  the  object  is  to  minimize 
the  trace  of  this  equation,  the  angles  0.^°’ 8 ^  and  n  could  be  con¬ 
sidered  as  control  variables  in  this  minimization.  Theoretical 
considerations  regarding  the  existence  of  a  minimum  with  respect 
to  .g^0’ and  n  are  made  in  Chapter  VI,  Only  experimental 
results  are  presented  here  which  might  give  some  insight  into  this 
optimization  problem.  The  sigma  for  the 
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IMUe  initial  misalignment  angles,  £^0,s!  listed  below  were 
s  o 

considered  first. 


0<?>s>  =(48.9,20,-22.5)*  (0,0,45)*  (0,45,0)*  (45,0, of 


To  make  a  comparison  a  =  was  plotted  for  the  four 

different  initial  IMUS  alignment  angles,  .  In  the  four 

cases  presented,  a  did  not  depend  on  the  initial  alignment  angles, 

—  o* S ^  indicatinS  that  a  minimum  of  the  sum  of  the  variances  of 

the  <r$f  does  not  exist  as  a  function  of  0^°’ s^.  It  was  clear  from 

the  results  that  the  individual  <r$.  do  depend  on  s),  so  that 

more  realistic  performance  criteria  might  be  to  minimize  the 

maximum  aA  as  a  function  of  ef?> s) .  Theoretical  difficulties 
vi  “° 

will  arise  with  such  a  criterion  because  the  maximum  principle 
cannot  be  used  to  obtain  necessary  conditions  for  optimality. 
However,  it  does  seem  that  the  equal  acceleration  orientation 
satisfies  the  above  criterion.  For  the  remainder  of  this  study, 
the  IMUg  orientation  is  such  that  O^  =  (48.9 , 20  ,  -22.5  )* 
(equal  acceleration  orientation)  was  chosen  because  the  three  o$. 
are  about  the  same,  and  there  is  no  reason  to  have  them  other¬ 
wise  at  this  point.  It  should  be  noted  that  the  random  accelera¬ 
tion  and  the  measurement  noise  were  assumed  to  be  such  that 
the  variances  on  each  component  were  equal,  i.e.,  each  com¬ 
ponent  of  velocity  difference.  In  the  case  in  which  these 
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variances  are  different  for  each  component  of  velocity  difference, 
the  initial  1MU  orientation,  0^o,s^  would  surely  affect  b . 

The  sigma  for  four  different  m  profiles  are  shown  in  Figure 
4.7.  The  four  different  nominal  trajectories  used  are  shown  in 
Figure  4.4,  and  are  essentially  such  that  p  approaches  its 
lower  value  of  55  deg  stepwise,  but  at  different  rates.  In  all 
cases  the  steady -state  cr£.  for  each  constant  portion  of  m 
(for  each  of  the  four  trajectories)  is  the  same.  Only  the  rate 
at  which  these  steady-state  values  of  o-^  are  approached  is 
different.  If  larger  and  were  used,  it  is  possible  that 
the  steady-state  values  would  not  be  reached,  and  that  the 
time  histories  for  the  four  nominal  trajectories  would  be  different. 
Actually,  the  size  of  the  step  determines  how  much  larger  the 
error  due  to  ^  in  the  velocity  differences,  y^0,s^,  is  than  that 
due  to  random  acceleration  and  measurement  noise.  The  larger 
this  difference,  the  more  confidence  the  estimator  has  in 
choosing  the  4^.  The  corresponding  minimum  variance 
estimates,  ,  for  two  cases  are  shown  in  Figure  4.8.  Further 
discussion  of  the  nominal  trajectories  are  given  in  the  next 

i 

j 


section 


SENSED  ACCELERATION  (FT/SEC) 


d)  Nominal  Trajectory  IB-6 


-3 

Note:  In  both  cases  A  =  .10,  °w  =  cr  =  ,  jj»  =  (.  5,  .  1,  -.3) 

Figure  4. 8b  Estimates  of  ^  Using  Nominals  IB-5  and  IB-6 


4.3  Experimental  Results  for  Initial  Misalignment  Angles 
and  Constant  Gyro  Drift  Rates  (Configuration  II) 


Effect  of  Random  Acceleration  and  Measurement  Noise 
on  Configuration  II 

Figures  4. 9  and  4. 10  show  the  effects  that  the  measurement 

noise,  cr  and  the  random  acceleration,  o\_,  have  on  the  esti- 
W  It 

mates  of  the  constant  gyro  drifts,  ci#  and  the  initial  misalign¬ 
ment  angles.  The  value  of  does  not  affect  the  *  too 

much  (Figure  4.9a),  except  decreases  faster  for  smaller  «r  . 

Xi  w 

The  steady -state  values  of  the  <r^  are  independent  of  o^,  as  is 

discussed  in  Section  4.5.  The  effect  of  on  the  is  stronger 

than  for  the  at  ,  especially  between  o_.  -  10"**  and  10 "4  ft/sec 
i  w 

(Figure  4. 9b) .  No  conclusions  about  the  steady-state  values  of 


still  decreasing  at  t  =  24  sec. 


The  effect  of  the  random  acceleration  («7p  )  on  the  estimation  of 

the  ^  and  is  illustrated  in  Fig’ires  4. 10a  and  4. 10b  respectively. 
-3  2 

For  <r_  =  10  ft/sec  ,  there  is  an  appreciable  steady-state  value 
K 

(~0.05  arc  min)  in  the  <r^,  that  is,  appreciable  when  compared  to 
the  error  due  to  the  constant  gyro  drifts,  As  shown  in  Figure 
4. 10b,  the  effect  of  the  random  acceleration  noise  levels  on  the 
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estimation  of  the  e.  is  considerably  more  pronounced.  For 

*r  =  10"5  and  10“6  ft/sec2,  the  ap  are  about  the  same  values. 

"i 

decreasing  significantly  to  about  0.0125  deg/hr  at  t  =  24.0  sec. 
These  results  indicate  that  it  would  be  desirable  to  keep  the 
measurement  noise  levels  such  that  1G"5  <  aw  s  10“4  ft/sec  and 
the  random  acceleration  noise  levels  such  that  10”5<  <rR  5  iq~4 
ft/sec^. 


Estimates  of  the  Initial  Misalignment  Angles 
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Figure  4-9b  Effects  of  Measurement  Noise,  o,^,  on  the 
Estimation  of  the  Constant  Gyro  Drift  Rates 


Figure  4-10a  Effect  of  the  Random  Acceleration,  cr^,  on  the 
Estimates  of  the  Initial  Misalignment  Angles 
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STANDARD  DEVIATIONS  IN  t  ST  (MATS  D  CONSTANT  GYRO  D#IFTS  (D€G/H») 


0  3  6  9  12  15  18  21  24 

TIME  (SEC) 

Figure  4. 10b  Effect  of  the  Random  Acceleration.  <rR,  on 
the  Estimates  of  the  Constant  Gyro  Drifts 

Effect  of  the  Nominal  Trajectory  on  Configuration  II 

The  trajectory  pitch  profile,  n ,  is  varied  to  study  its  effect  on 
the  estimation.  The  variations  are  such  that  the  effects  of  the 
rate  of  change  of  u,  total  amount  of  change  in  n,  and  direction 
of  the  change  in  n  on  the  estimates  can  be  observed.  It  should 
be  noted  that  the  constant  drift  rates,  do  not  appieciably 
affect  the  estimates  of  the  initial  misalignment  angles,  ^ , 
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because  the  errors  due  to  ei  are  much  smaller  than  those  due  to 
the  *1  for  short  estimation  times.  For  this  reason  the  conclu¬ 
sions  regarding  the  ^  deduced  from  this  section  can  be  used  in 
the  previous  section. 

Figure  4. 11  shows  six  trajectory  pitch  profiles  in  which  they 

change  as  ramp  functions.  The  corresponding  standard  deviations 

in  the  4^,  and  the  standard  deviations  in  the  <x^,  are 

shown  in  Figures  4. 12a  and  4. 12b  respectively.  It  is  seen  from 

Figures  4. 12a  and  4. 12b  that  the  higher  rate  of  changes  in  TP 

result  in  faster  decreasing  <r&  and  <ta  .  However,  the  steady 

i  ei 

values  in  the  <t6  ,  that  is  at  c  =  24.0  sec,  are  about  the  same 

(about  0.025  arc  min).  The  <r|  are  still  rapidly  decreasing  at 

i 

t  =  24. 0  and  the  rate  of  decrease  of  the  is  proportional  to  the 

rate  of  change  of  the  pitch  angle.  For  TP-2f,  the  oa  ,  have 

i 

decreased  from  0. 15  deg/hr  to  about  0.040  deg/hr. 
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Figure  441  Trajectory  Pitch  Angle  Profiles  for  Nominal 
Trajectories  II-2a  through  II-2f 


Figure  442a Standard  Deviations  in  the  Estimated  Misalignments,^ 
for  Nominal  Trajectories  II- 2a  through  II-2f  V 


STANDARD  OrVIArtONS  IN  ESTIMATED  CONSTANT  GYRO  DRIFTS  (DEG/HR) 


Figure  4d2b  Standard  Deviations  in  the  Estimated  Constant  Gyro 
Drifts,  o-A  ,  for  Nominal  Trajectories  II-2a  to  U-2f 


Next,  the  effect  of  the  total  change  in  TP  on  the  estimates  is 

considered.  The  pitch  profile  is  changed  discretely  at  t  =  10  sec, 

as  shown  in  Figure  4, 13at  These  nominal  trajectories  are 

referred  to  as  Cases  H-3a  to  II-3g.  The  corresponding  tr£  and 

*i 

are  shown  in  Figures  4. 13b  and  4.13c,  d  respectively.  The 
rate  of  decrease  of  the  sigmas  for  ^  and  ^  for  the  various 
"step”  pitch  profiles  is  proportional  to  the  larger  changes  in  TP. 
This  is  reasonable  because  the  larger  the  change  in  the  pitch 
profile,  the  larger  the  contribution  of  the  and  to  the 
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velocity  differences,  v^0>  s^,  as  compared  to  the  contribution  due 

to  the  random  acceleration  and  measurement  noises.  From 

Figure  4, 14a  it  is  noticed  that  <r£  and  cr*  improve  when  TP 

12“ 

increases  (!I-3e  and  H-3g)  and  it  is  relatively  insensitive  to  down¬ 
ward  changes.  On  the  other  nand,  Figure  4. 14b  shows  that  cr£ 

3 

improves  as  TP  decreases  (II-3d  and  n-3f)  and  does  not  change 
much  when  TP  increases  (II-3g  and  n-3e).  This  indicates  that 
a  combination  of  an  upward  and  a  downward  maneuver  would  be 
better  for  extracting  the  constant  gyro  drift  .  This  is  further  in¬ 
dicated  by  the  results  shown  in  figure  4. 15b  using  nominal  li-5c. 


mi  fisc) 

Figun 4.13a  Nominal  Trajectories  II-3a  through  II-3g 
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STANDARD  DEVIATIONS  IN  ESTIMATED  CONSTANT  GYRO  DRIFTS 
and  (DEG/hR) 


Figure  4. 13b  cr^  for  Nominals  II-3a  through  H-3g 
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Figure  4. 13c  <t*  and  ar£  for  Nominals  n-3a  through  II-3g 
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STANDARD  DEVIATIONS  IN  CONSTANT  GYRO 
DRIFT,  of  (DEG/HR) 


Figure  4. 13d  Standard  Deviations  in  the  Estimated  Constant  Gyro 
Drifts,  ,  for  Nominals  H-3a  Through  n-3g 


The  trajectory  is  next  varied  in  such  a  way  that  the  step  of  10  deg 
(trajectory  TP  H-3f)  is  approached  at  different  rates.  This  type  of 
variation  gives  the  tradeoff  in  performance  of  Configuration  n  as  a 
function  as  less  abrupt  changes  in  the  pitch  angle.  The  trajectory 
pitch  angles  are  shown  in  Figure  4.14a.  The  corresponding  0$  and 
are  shown  in  Figures  4.14b,  4.14c  and  4.14d,  respectively.  As 
would  be  expected,  the  trajectory  profile  with  the  most  abrupt 


changes  produces  better  results,  for  as  the  cr^  steady  out 
within  one  second  and  the  most  abrupt  changes  in  the  trajectory  con¬ 
tribute  a  larger  amount  to  v^°’ 8  ^  for  a  longer  period  of  time.  For 
the  the  slower  ramps  produce  better  results  because  the  <r^ 
take  longer  to  steady  out,  and  better  results  would  be  obtained  if 
the  changes  are  made  later. 
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It  should  be  noted  that  if  the  estimation  time  is  long  enough,  it 
doesn't  seem  to  matter  how  the  trajectory  changes  -  as  long  as 
it  does  change.  This  is  true  for  the  at  least,  because  the 
all  seem  to  steady  out  to  the  same  values.  This  could  be 
explained  in  terms  of  c rw  and  <r  and  the  discussion  of  Section  4.5 
Note  that  the  slower  the  ramp  decreases,  the  better  the  ?  are 
estimated  (contrary  to  the  situation  for  the  * ).  This  is  clear 
from  Figure  4.14c  and  4.14a.  Also,  is  better  than  £  and  &n 
because  the  trajectory  pitch  angle  is  decreased  from  65  deg. 
This  behavior  was  observed  for  the  step  profile  of  Nominal 
Trajectories  H-3a  to  H-3g. 

Figures  4. 13  and  4. 14  indicate  that  <xa  is  better  than  a&  and  aA 

€1  k2  €3 

if  m  changes  downwardly  from  65  deg,  and  vice-versa.  Figure 
4. 15b  shows  the  for  a  combination  of  a  downward  and  an 
upward  maneuver.  The  specific  trajectories  used  are  C-5a, 
b,  and  c  shown  in  Figure  4. 15a  with  cr^  "  =  10"*.  It  is 

clear  from  Figure  4.15b  that  the  can  be  equally  well  estimated, 
and  is  less  than  0. 055  deg/hr  after  24  sec  of  estimation. 

If  the  trajectory  pitch  profile  were  varied  in  an  optimal  manner, 
better  estimates  would  surely  result. 
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Figure  4. 14a  Nominal  Trajectories  H-4a  Through  II-4e 
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Figure  4. 14b  Standard  Deviations  in  the  Estimated  Misalign 
ment  Angles,  crA ,  for  the  Various  Ramps  of 
Nominal  Trajectories  H-4a  Through  D-4e 
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Figure  4.14c  Standard  Deviations  in  the  Estimated  Constant  Gyro 
Drifts,  and  ,  for  the  Step  Profiles 
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Figure  4. 14d  Standard  Deviations  in  the  Estimated  Constant  Gyro 

Drifts,  cr*  ,  for  the  Step  Profiles  (Nominal  Trajectories 
H-4e) 


4.4  Experimental  Results  for  Initial  Misalignment  Angles 
and  Mass  Unbalance  Drift  Rates  (Configuration  HI) 

The  effect  of  the  nominal  trajectory  and  the  initial  orientation  of 

IMUS  on  the  estimates  has  been  discussed  in  Sections  4.2  and  4.3. 

These  effects  on  Configuration  m  are  not  studied  here,  but  it  is 

expected  that  the  results  obtained  in  the  previous  sections  would 

generally  apply,  at  least  for  the  estimation  of  The  estimation 

of  the  kc  and  kT  errors  would  possibly  require  more  trajectory 
bi  Ai 

maneuvering. 

The  first  set  of  standard  deviation  curves  (Figures  4. 17  and  4. 19  ) 

show  the  effects  of  the  random  acceleration  on  the  estimates.  The 

sigma  4.  are  shown  in  Figure  4. 17,  and  it  i*-  clear  from  these 

curves  that  the  trajectory  pitch  angle  could  have  been  changed 

earlier,  since  the  steady-state  values  of  about  0.67  are  reached 

rather  quickly.  Estimates  of  the  ^  such  that  <4^  <  0. 1  arc  min 

are  obtained  in  all  cases,  except  for  crR  =  10  ft/sec  .  The 

spin-axis  mass -unbalance  standard  deviations,  <4  ,  for  the 

different  cr_  are  shown  in  Figure  4. 19a.  Except  for  the  case 
JK 

where  a  =  10"^,  p  could  also  have  been  changed  earlier.  The 
input-axis  mass-imbalance  standard  deviations,  o-j^  ,  are  shovvn 
in  Figure  4. 19.  These  estimates  are  somewhat  better  than 
those  for  the  spin-axis  terms.  Again  the  pitch  angle  could  have 
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been  changed  earlier,  say  at  t  =  5.0  sec  .  These  results  show 

that  if  lCf4  ft/sec  and  cr^  *  10“4  ft/sec  (except  the  kg^  terras 

need  <  10'4)  and  estimating  times  aie  in  excess  of  20  sec, 
ft 

excellent  estimates  of  4^,  kg^,  and  kj^  would  be  obtained.  Based 
on  these  results  it  is  conjectured  that  a  well  chosen  maneuver 
would  improve  the  estimates,  especially  the  estimates  of  the  spin- 
axis  terms,  kg. . 

Next  the  effect  of  the  measurement  noise  on  the  estimates  is  con¬ 
sidered.  The  nominal  trajectory  is  TP-IH-3  in  Figure  4. 16. 

This  is  essentially  the  same  as  Nominal  Trajectory  m-2,  except 
that  TP  is  changed  from  6*5  deg  at  5  sec  instead  of  at  10  sec.  The 
corresponding  estimation  sigmas  are  shown  in  Figures  4.  IS  and 

C 

4.20.  Typical  minimum  variance  estimates  (for  -  10  and 
°R  "  10”4)  are  shown  in  Figure  4.21.  Generally,  the  measure¬ 
ment  noise  does  not  affect  the  estimates  as  in  the  case  of  the 

random  acceleration  a^,.  Figures  4. 18  and  4.20  show  that  there 

-3  -4 

is  not  much  difference  in  the  sigmas  for  cr^  =  10“  and  <rw  =  10  . 


Figure  4.16  Nominal  Trajectories  331-2  and  TP  HI- 3 
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STANDARD  DEVIATION  IN  THE  MASS  UN&/  1.ANCE  DRIFT  RATES,  i5  (OEG/HR/G)  >T)  STANDARD  DEVIATIONS  IN  THE  ESTIMATED  MISALIGNMENT  (ARC  MIN) 


igure  4.17  Effect  of  Random  .•  ^deration  on  the  Estimates 
of  Initial  Misalignments,  dr. 


Figure  4.iS  Effect  of  Measurement  Noise  on  the  Estimation 
of  the  Initial  Misalignment  Angles,  i}>. 


STANDARD  DEVIATION  IN  THE  MASS  UNBALANCE  ORI 
RATES,  £,  (DEG/HR/G) 


Figure  4. 20b  Effect  of  on  the  Estimation  of  kj 
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Figure  4,21  Minimum  Variance  Estimates  of  and  kj„ 


i=l,2, 3. 


4. 5  Discus 3 ion  of  the  Steady  State  Estimates 

The  knowledge  of  the  steady  state  estimates  and  of  the  error  bounds 
on  the  estimation  covariances  is  not  required  directly  in  the  maneu¬ 
ver  optimization  problem.  However,  this  type  of  discussion  could 
help  to  answer  certain  problems  of  controllability  which  might  arize. 
For  example,  the  final  error  covariances  which  need  to  be  speci¬ 
fied  for  certain  trajectory  optimization  problems  (Section  5.4) 
could  not  be  less  than  the  steady  state  results  given  below.  For 
this  reason  some  results  which  might  be  applied  to  the  asymptotic 
estimator  of  the  covariances  are  included  [  4  ].  To  present  the 
main  results  of  this  section,  the  general  state  and  measurement 
equations  are  written 

=  0k,  k-1  Be,  k-1  *  2k- 1 

Be  =  M*k  +  Ik  • 

Cov  (w)  =  Q  and  Cov  (v)  h  R 

We  write  M  instead  of  Mk  in  the  last  equation  because  M  =  ( I  :  O), 
a  constant  matrix  in  the  present  problem. 

Definition:  The  above  system  is  said  to  be  q-stage  observable 
(1  s  q  s  N)  on  an  interval  *o  *  «k  "  fcN  if  and  only  if  the  matrix 
Mkjk-q+  1  is  defined  below)  is  positive  definite  for  arbitrary 

tfc  and  q  such  that  tj  s  tk„q+  i  and  tk  «  tN. 
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m*V!mv 


In  the  present  problem 


;2*5Ilo2*«2 


which  can  be  made  positive  definite  for  the  three  configurations 
considered  here. 


Result  1:  If  the  process  noise  is  zero  and  if  the  system  is  q-stage 
observable,  then  the  error  covariance  matrix  of  the  estimates  Pk 
vanishes  as  k  if  O  more  rapidly  than  l|$k  ^ 

increases,  where 


Result  2:  If  the  process  noise  is  zero  and  the  system  is  q-stage 
observable,  then  Pk  becomes  essentially  independent  (that  is,  for 
k  large  enough)  of  the  a  priori  statistics  P0 . 


Result  3:  If  the  measurement  noise  is  zero,  then  Pk  satisfies 

(a)  MPjj  =  O  and  PkMA  =  O  for  each  k, 

(b)  Pk  is  nonnegative  definite,  but  never  positive 
definite 

(c)  If  Pk=  k*j  Pk-1  *k  k->1  +  Qk.1  is  positive  definite, 
and  m<n  (where  M  is  mxn),  then  P^can  never  vanish. 
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Further,  the  error  in  the  estimate  x^  -  x^  -  is  any  element  of 
the  null  space  of  M,  i,  e. ,  -  q  for  arbitrary  k. 

Result  4:  The  error  bounds  for  Pjj  for  the  q-stage  observable 
system  are  given  by 


where  p£  is  any  nonn'igative-definite  matrix  (never  positive  defi¬ 
nite  and  generally  nonzero)  and 

Wk,l  ^  *k,iQi-l  %,i 

It  would  be  useful  to  compute  the  above  error  bounds  as  a  function 
of  the  various  system  parameters. 


5 


§ 


If  it  is  assumed  that  estimator  IB  is  q-stage  observable,  then  the 

i 

above  results  may  be  used  to  interpret  Figures  4.4  and  4. 5.  Figure 
4.4  illustrates  that  all  of  the  curves  seem  to  be  going  to  the  same 
steady-state  values  for  Vp  fixed  and  varying  <rw .  This  observation 
is  implied  by  Result  1  above.  Result  2  implies  that  the  estimation 
covariance  will  become  essentially  independent  of  the  initial  statis¬ 
tic  P  although  no  experimental  results  are  presented  as  verification. 
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Figure  4. 5  indicates  that  the  <ta  will  reach  steady-state  values 
which  are  proportional  to  crR.  This  is  implied  from  Result  3(c) 
above.  Figure  4. 4  indicates  that  the  estimates  are  not  improved 
for  decreasing  aw  if  the  estimation  time  is  sufficiently  long.  These 
results  are  not  used  to  discuss  Figures  4.9,  4. 10  and  Figures  4. 17, 
4„  18,  4. 19,  4. 20,  because  steady-state  values  were  not  reached 
for  these  configurations.  It  should  be  noted  that  the  bounds  on 
the  covariance  equation  which  are  given  in  f  4  ]  for  the  discrete  case 
are  similar  to  those  given  in  Section  15  and  16  of  f  5  ]  for  the  con¬ 
tinuous  case. 

Remark  4  ♦  1  The  experimental  results  of  this  chapter  indicate 
that  significant  improvements  in  the  identification  of  the  specific 
error  parameters  considered  can  be  obtained  by  changing  the 
nominal  trajectory.  Furthermore,  certain  changes  result  in 
estimation  variances  which  decrease  quicker  than  for  other  changes. 
The  best  time  to  make  a  trajectory  change  for  a  particular  error 
parameter  seems  to  be  when  the  cr  for  this  error  is  approaching  a 
steady  state  value.  Section  4. 5  indicates  that  if  the  system  satis¬ 
fies  certain  properties  the  steady  values  (that  is,  c— *«)  can  be 
estimated  independently  of  the  nominal  trajectory.  In  section  5. 1 
it  will  be  seen  that  some  of  these  properties  are  related  to  the 


"no-noise"  situation 


CHAPTER  V 


ANALYSIS  OF  TRAJECTORY  MANEUVER  AND  FORMULATION 
OF  THE  OPTIMAL  CONTROL  PROBLEM 


5. 1  Identification  of  errors  under  ideal  conditions 

For  the  purposes  of  determining  the  minimum  number  of  trajectory 
chants  for  the  uniqueness  of  the  solutions,  we  consider  the  mini¬ 
mum  number  of  trajectory  maneuvers  which  are  required  to 
identify  the  IMU  error  parameters  under  ideal  conditions,  that  is, 
under  the  conditions  of  no  random  disturbances  (°r  “  °w  ~  °)‘  To 
make  the  discussion  as  uncomplicated  as  possible,  the  simplest 
configuration  (Configuration  IB )  is  considered  first.  For  the 
no-noise  situation  we  have 

jv^(k+l,k)j*  =  k*1  ^v(°>s)(k+i)  -  v^°,s^(k  )j  (5.1a) 

If  all  subscripts  are  dropped  and  k  takes  on  two  different  values, 
the  above  equation  gives 

=  Vf  and  A g  i  =  v  2  (5.1b) 

where  the  subscripts  1  and  2  correspond  to  kj+  1  and  k2+  1 
respectively. 
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Property  5. 1  For  the  no-noise  situation,  at  least  one  trajec¬ 
tory  maneuver  must  be  made  in  order  that  the  misalignment  angles 
be  determined  uniquely. 


Demonstration  The  matrix  Aj  is  of  the  skew  symmetric  form 


Ai  = 


0 

a3 

-a2 


-a3  a2 

0  -d-i 

H  0/ 


for  which  the  eigenvalues  are  =  0 ,  Xg  -  +j  (fall  and  Xg  =  -j  Itafl , 
whpre  !la||  is  the  magnitude  of  the  sensed  acceleration.  This 
means  the  three  dimensional  space  in  which  *  is  defined,  £3 ,  can 
be  spanned  by  the  set|x* ,  x^ ,  x^  where  Xj  is  the  eigenvector 
corresponding  to  the  eigenvalue  x* ,  and  any  vector  x  in  Eg  can  be 
represented  by 


12  3  19  —2 

2  =  +  “gX  =  OjX  +U^x“+^x 


so  that 

Ail  =  Xx  q^x*  +  Xg  Q^x^+  Xg  ~  *2  e2-^+  *2 

(5. Id) 

If  *  has  a  component  in  the  null  space  of  x^  =  0 ,N0  (i.e. ,  NQ= 

{x  in  E3  1  Ax  =  0}  ),  then  ±  cannot  be  solved  for.  On  the  other 
hand,  if  ±  does  not  have  a  component  in  NQ ,  then  any  multiple  of 
a  vector  in  N0  can  be  added  to  ±  and  this  will  also  satisfy  the 


equation  =  v. 


For  convenience,  the  common  factor  which  takes  care  of  the  thrust 
acceleration  can  be  factored  from  the  A  matrix,  so  that  we  have 


Ai*  =  Zi 


(5,2a) 


where 

0  -t3  t2 

t3  0  -tj 
t2  ti  0 

The  eigenvalues  of  A^  and  A2  are  thrn  the  same,  namely, 
=  0,  Xg  =  vQ  and  Xg  = 


If  the  estimation  were  performed  at  the  same  time,  then  the 
problem  is  the  following:  find  the  vector  w  such  that 


(5.2c) 

has  rank  3,  then  a  unique  ±  may  be 


obtained  by  multiplying  Equation  (5.2c)  and  then  computing  ± 


by  the  equation 


(5. 2d) 
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When  the  estimation  is  performed  recursively,  Aj$  =  yj  must  be 
solved  for  %  first,  and  then  A2±  =  y£  is  solved.  The  vector  ± 
may  be  represented  by  the  eigenvectors 


±  =  x  j+  ofg  2^1  =  a\  x{  +  a2  (5'3a) 


The  vector  ^  is  then  given  by 


h 


5)- 


(5.3b) 


1  2 
The  coefficients  01 2  may  b®  solved  for  in  terms  of  y^  and  Xj,  which 

are  both  known.  This  corresponds  to  the  solution  of  minimum  norm 
(pseudo-inverse).  If  this  solution  is  denoted  by  iq,  then 
i  -•  ^  where  Cj  is  ar.  unknown  constant  and  xj|  is  the  eigen¬ 

vector  oi  Aj  corresponding  to  the  zero  eigenvalue.  Next  the 
Equation  A2*  =  £2  is  solved  for  a  corresponding  *2  >  311(1  £  is  £herl 
given  by  4  =  ^2  +  c2xf .  These  Equations  can  be  solved  for  q  and 
C2,  and  thus  give  ±  uniquely. 


Remark  5. 1  When  random  aistrubances  are  added,  a 

similar  discussion  goes  through  with  the  reasoning  that  the  change 
inu(and  therefore  the  change  in  the  matr  ix  A)  must  be  large  aiough 
so  that  the  error  due  to  ±  can  be  detected  over  the  errors  due  to 
the  random  disturbances. 
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Remark  5,2  The  problem  of  uniqueness  of  the  estimates 

for  Configuration  II  can  be  considered  as  in  the  above  discussion. 
In  this  case  the  matrices  A$  are  3x6  instead  of  3x3,  and  it  is  re¬ 
quired  that  [A*  1  —  |  A*]*  have  rank  6.  The  matrix  A±  is 

the  3x6  matrix  corresponding  to  the  i*h  observation,  and  is  of  the 
form 


*0 

-t3 

*2 

1 

i 

1 

"0 

-fc3 

fc2 

fc3 

0 

-*1 

*  P 
1 

fc3 

0 

-tl 

:fc2 

h 

0. 

l 

1 

:fc2 

tl 

0. 

where  a,  p.  and  the  tj  change  with  time.  With  r.o  trajectory 
maneuvers,  the  matrix  A*A,  where  A*  =(A*  j  A^]  has  two  eigen¬ 
values  which  are  essentially  zero  and  one  which  is  close  to  zero, 
even  though  the  factors  a ,  p  change  at  different  rates.  With  one 
maneuver,  the  matrix  A* A  still  has  the  same  zero  eigenvalue 
properties,  and  with  two  trajectory  maneuvers,  the  matrix  A*A, 
where  A  =  [Aj  j  Ag  j  A3] ,  has  no  zero  eigenvalues.  Thus  it  is 
sufficient  for  configuration  II  to  have  two  trajectory  maneuvers  in 
order  that  and  t  be  identifiable  under  ideal  conditions. 

Remark  5.3  The  uniqueness  property  for  Configuration  HI 

is  again  discussed  as  above,  where  the  matrix  Aj  is  now  of  the 


form 


I 


'o 

“*3 

*2 

i 

i 

■»» 

0 

“*32  *23 

0 

‘*33 

*21 

*3 

0 

-fei 

1  8 
i 

*31 

0 

1 

“*13  ' 

t 

*32 

0 

-*11 

r*2 

tl 

o_ 

* 

i 

:*21 

*12 

0  • 

1 

”*22 

*13 

0  . 

where  a,  #,  t^,  ty  ,  change  with  time.  In  this  case  it  is  required 
that  [Aj  j —  jA*]*  be  of  rank  9.  For  this  configuration  A  A  has 
zero  eigenvalues  when  no  maneuver  is  executed,  and  no  zero  eigen¬ 
values  when  three  maneuvers  are  executed. 

From  Property  5. 1  and  the  above  remarks,  the  following  conjecture 
is  made 

Conjecture  5. 1  In  order  that  the  parameters  of  a  particular 

parameter  identification  problem  be  identifiable  under  ideal  con¬ 
ditions  it  is  sufficient  that  a  minimum  number  of  trajectory  man¬ 
euvers  be  executed.  The  number  of  trajectory  variations  can 
be  obtained  from  the  transition  matrix  of  the  system.  In  the 
presence  of  random  disturbances  It  is  plausible  that  a  similar 
situation  prevails,  with  the  maneuvers  being  large  enough  to  detect 
the  changes  in  the  errors  due  to  the  error  parameters  over  that 
d?ie  to  the  random  disturbances. 

Remark  5.4  The  lower  bounds  of  the  minimum  number  of 

trajectory  maneuvers  could  be  useful  in  determining  the  ’’optimal 
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controls*'  discussed  below  and  in  Chapters  VI  and  VII.  In  Chapter 
VI  it  is  shown  that  the  optimal  controls  are  often  bang-bang.  A 
lower  limit  on  the  number  of  switches  in  the  controls  reduces  the 
search  procedure  which  might  be  used  to  obtain  good  initial  guesses 
for  the  computation  of  the  optimal  controls  as  discussed  in  Chapter 

vn. 

Figures  5. 1  and  5.2  illustrate  the  standard  deviations  in  the  esti¬ 
mated  misalignment  angles  and  the  corresponding  estimates,  ^ , 
respectively,  for  =  <rw  =  0.  The  trajectory  pitch  angle  profiles, 
ju,  for  these  two  cases  w-ere  such  that  u  =  constant  (IB-1)  and  p 
changed  discretely  downward  by  5  degrees  at  0.2  seconds  (IB-2). 
The  "squares"  in  Figure 5.  lb  indicate  the  estimates  for  the  second 
trajectory,  where  the  actual  values  for  ±  are  (0.5,  0.1,  -0.3). 

The  variances  and  the  estimates  for  Configuration  II  are  shown  in 
Figures  5.2a  through  5.2b.  With  no  process  and  measurement 
noise,  the  standard  deviations  in  the  estimated  misalignments, 

0“$. ,  do  not  drop  immediately  to  zero,  but  seem  to  steady  out 
at  approximately  0. 67.  At  t  =  0. 4  sec  the  pitch  profile  changes 
discretely  by  5  deg,  and  these  rigmas  then  go  essentially  to  zero. 
The  estimated  misalignment  angles,  are  estimated  perfectly 
after  t  =  0.4  sec.  The  reason  that  these  angles  are  not  given  cor¬ 
rectly  before  t  =  0.4  sec  is  that  the  estimated  misalignments  need 
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TIME  i$EO 


Figure  5.  la  Initial  Misalignment  Sigmas,  cr§  ,  for  the 
Nc-Noise  Situation  orR  -  aw  =  0  i 


<  -0.20 


TIME  (SEC) 


Figure  5.  lb  Estimates  of  the  Misalignment  angles,  for 

the  No-Noise  Situation  =  °\y  =  0  factual  £  values 
are  0.5,  0.1,  and -0.3). 


TIME  (SEC) 


Figure  5. 2d 


Estimates  for  the  case  =  ffw  =  0 
(Actual  Values  were  0, 15,  0. 05  and  -0. 15) 
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not  be  unique.  With  the  estimation  of  the  constant  drift  rates, 
there  is  some  difficulty  at  first.  At  0.7  sec  it  seems  that  there 
is  a  "numerical  instability"  in  the  However,  the  sigmas 
settle  down  to  essentially  zero  at  t  =  1.0,  and  the  estimated  ej 
come  out  perfectly. 

Using  the  nominal  trajectories  shown  in  Figure  5.3  below  (Nominals 
m-la  and  HI- lb),  the  sigmas  for  the  misalignment  angles,  <7$,,  and 
the  sigmas  for  the  mass-un.balance  gyro  drift  rates,  and  <7^, 
far  the  no-noise  conditions  ( =  0)  are  shown  in  Figures  5.4 
and  5.5  respectively.  The  corresponding  estimates  and 

,  are  shown  in  Figure  5.6.  It  is  noted  in  Figures  5. 4  and  5. 5 
that  the  sigmas  do  not  go  to  zero  for  the  constant  pitch  profile 
(TP  m-ia).  When  the  pitch  profile  is  varied  (TP  IH-lb),  the 
sigmas  go  down  to  zero.  For  this  particular  case  perfect  esti¬ 
mates  of  are  obtained  for  both  nominal  trajectories  (Figure  5.6  ). 
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ft  (ARC  MIN) 


0  0.2  C.4  0.6  0.8  i.O 


TIME  (SEC) 

Figure  5.4  Misalignment  Sigmas,  o-^.,  Using  Nominals 


a.  SPIN  AXIS  DRIFT  RATE  SIGMAS 


0  0.2  O.A  0.6  0.6  1.0 


TIME  (SEC) 

b.  INPUT  AXIS  DRIFT  RATc  SIGMAS 

Figure  5. 5  Mass -Unbalance  Drift  Rate  Sigmas, 

ando^  using  Nominals  IH-la  and  Ill-lb  with 
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a.  MISALIGNMENT  ANGLE  ESTIMATES,  £  (ACTUAL  VALUES  0.3,  0,  -0.3) 


b.  SPIN  AXIS  DRIFT  RATE  ESTIMATES,  K$  (ACTUAL  VALUES  0.15,  0  JO,  and  -0.05) 


0  0.2  0.4  0.6  0.8  1.0 

TIME  (SEC) 


5.2  Formulation  of  the  Optimal  Control  Problem 


In  this  section  the  manner  in  which  the  control  (trajectory  pitch 
angle)  emers  the  covariance  equation  for  the  error  parameters  « 
is  indicated,  and  the  pertinent  optimization  equations  are  formu¬ 
lated.  Under  the  assumptions  for  Configurations  IB,  n  and  HI, 
the  state  and  observation  equations  are,  respectively 
x  =  Ax  +  Bw 


z  =  Mx  +  v 

where  x  =  (v^0»s),«)*,  a  (3  +  N)  x  1  vector,  where  N  indicates  the 
number  of  components  in  a ,  =  *  (Configuration  I),  «jj  =  (i,je )* 

(Configuration  H),  =  (!>  kg’-I^*  (Configuration  ID),  w  is  a  3x  1 

vector  of  random  (uncorrelated)  vehicle  accelerations,  and  v  is  a 
3x1  vector  of  uncorrelated  observation  disturbances.  The 
respective  covariances  of  these  random  vectors  are 

Q  =  ccv(w)  =  q  -  I  and  R  =  cov(v)  =  r  *  I 


The  more  convenient  notation  q  and  r  are  used  from  here  on  instead 
of  and  o-^..  The  matrix  B  is  defined  by  B  =  ^T-°-S-)  J  ?-Bl’ 
^(3+N)  x  3^,  where  T^°5  is  the  3x3  coordinate  transformation 
from  the  master  to  the  slaved  IMU.  This  transformation  involves 
the  angles  which  orient  the  slaved  IMU,  8^°’ s),  i  -  1, 2, 3,  which 


no 


are  now  assumed  to  be  fixed  for  the  "equal  acceleration”  orienta¬ 
tion  and  T^m,°^  is  the  3x3  coordinate  transformation  from  the 
vehicle  to  the  master  IMU  coordinates.  This  transformation 
involves  the  trajectory  pitch  profile,  0^’ which  we  write  as 

9(2°’m)(t)a(l0(t)+  „(t) 


to  denote  a  nominal  profile  vQ(t)  and  an  off-nominal  control  in  the 
profile  p(t>.  The  observation  matrix  M  is  given  by 


M 


=  I  !  O  (3x 

.  I 


(3  +  N)> 


The  matrix  A  is  given  by 


A  = 


O  !  Ai 
O  '  0 


i  =  i.n,ra, 


where 


Aj  -  k*[aiS)] 
An  =  [k*[4s)] 
Am=  [k*[4s)] 

It  is  recalled  that  the  matrix 


(3x3) 

;kd.r.  t[4S)j] 

(3x6) 

ikm.u.[4S)]^[4S)]2dT] 

(3x9) 

£a^s|J  is  defined  by 

ill 


/ 


[*>]- 


V 


o 

4 


-4S> 


3 

(s) 

2 


0 


(s) 

where  ag.,  i=l,2,3  are  the  components  of  the  nominal  acceleration 

(g\ 

in  the  IMUS  coordinates,  1  •  The  superscript  and  subscript  s 
on  the  acceleration  will  be  dropped,  unless  the  meaning  is  not  clear. 
In  particular,  each  component  is  given  by 
ai(t)  =  TjU)  a{“>(t) 

where 


T  a 


T., 


LT3j 


=  T(°’s)  T(m’0)  ‘  \  0 


For  the  equal  acceleration  orientation  and  0^,m^  =  -25°, 

T  h  =  -^(1, 3,  1)*.  The  thrust  acceleration  a^m^(t)  is  of  the 
form  a/(P  -t),  which  approximates  a  vehicle  with  a  constant  pro* 
pellant  burning  rate. 

The  matrix  j  is  defined  by 


[■ 


ks 


>2 


"al 

0 

0 

i  H 

0 

0" 

0 

a2 

0 

!  o 

•3 

0 

_0 

C 

a3 

!  c 

C 

al_ 
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when  multiplied  out, 


dt  = 


0  *a32  **23  J  0  '**33  *21 

Hi  0  *a13  !  *32  0  'all 

"a21  a12  ®  i  “a22  a13  ® 


where 


s  Tj(t)  a(8m,(t> 


fV>4m) 

n 


(r)dr,  1,  j  «  1, 2,3. 


If  we  make  a  small  angle  approximation  for  /*{t),  we  obtain 


*$»> ‘  FI 


fCOoSi 

\-sin  j^/ 


-sin 
+  u(t)  [  0 

cos^> 


I 


For  ^(t)  =  -25°  and  the  equal  acceleration  orientation. 


af s\t )  =  -2 L. 
W  0_t 


0/  .81 5\ 

+  "(t>  H 

\  -.409/ 


Thus  the  control  p(t)  enters  linearly  iiito  Aj,  An,  and  in  a  compli¬ 
cated  bilinear  fashion  in  Ajjj.  That  is, 

3ij(t)  =  (ci(t)  +  d^tJ^lt^J1 (Cj(t)  +  dj(r)  #i(r)  dr^ 


Remark  It  is  noted  that  vehicles  with  variable  thrust  engines  (for 

example,  airplanes  and  ships)  can  be  considered  in  this  framework 

as  well.  In  this  case,  the  term  a^m^  becomes  a  control  variable, 
(s) 

and  jig  is  of  the  form 


bujxj* 


i 

t 


^  =  ux(t)  cx(t)  + c2{t)  ux(t)  U2(t), 
where  cx(t)f  c2(t)  are  known  3x1  vectors,  ux(t)  is  the  thrust  ac 
celeration,  and  U2(t)  is  the  trajectory  pitch  angle. 


The  estimation  covariance  equation,  P  =  E  [(x  -  x)(x  -  x)*] ,  is 
specified  by 

P  =  AP  +  PA*  -  PM*FT*MP  +  Q 


If  P  is  partitioned  in  the  obvious  way 
P  = 


pXX  pX® 


potx  p 


aa 


there  results 

pxx  =  pox*  A*  +  AjI^-r^P^P^+q-  I 

pax  =  ^  Aj  -  r"1  P"x  P*^ 
pF1*  -  _r-f  pPfX  p»x* 


Remark  If  the  column  vectors  of  PQX  are  written  Px » •  •  •  *  PN>  411611 


/p*Pl...pfPN\ 

j/JrXpflfX  _  :  ; 

V^v'S  w 


and 


Trace  P 


<*'?>*-£  pfPr  f  i(rf)2 


=1  i= 1  j=l 


Thus 


Trace  P°°U)  =  £  Pj(°)~^  fNpjjrfdt 

i=i  i=iFr° 
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so  that  minimizing  Trace  P^t)  is  the  3ame  as  maximizing 


The  maximization  is  performed  with  respect  to  changes  in  the 
trajectory  pitch  angle,  p(t). 


5.3  Trajectory  Constraints 

The  constraints  on  the  vehicle’s  trajectory  will  determine  con¬ 
straints  on  the  pitch  angle  p.  In  the  simplest  forms,  the  trajectory 
constraints  might  be  specified  in  several  different  ways  as  listed 
below. 


!h(t)| 

<  Mj 

(C-l) 

<  m2 

(C-2) 

j£%(t)dt<M3 

(C-3) 

li  “<t)d 

t|5M4 

(C-4) 

and  various  combinations  of  the  above  constraints.  The  constraint 
(C-4)  is  the  most  difficult  to  treat  analytically*  and  will  not  be 
discussed  further  here.  Constraint  (C-2)  can  be  considered  by 
defining  A  as  the  control,  and  u  as  a  state  variable.  If  constraint 
(C-l)  is  also  in  force,  then  the  problem  becu.aes  one  with  bounded 
state  variables.  In  this  case  the  specification  of  the  necessary 
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ccno  'Jons  which  the  optimal  control  £  must  satisfy  become  quite 
complicated.  However,  with  some  of  the  cost  functionals  in  the 
next  section  it  can  be  shown  (Chapter  VI)  that  with  (C-l)  in  force, 
we  require  }p(t)|  =  and  with  (C-2)  in  force  we  require  !£(t)!  =  Mg, 
that  is,  n  and  m  are  bang-bang  controls.  Thus  with  (C-l)  and 
(C-2)  in  force,  then  it  is  plausible  thatM(t)  =  until  ^(t)  =±M| 
and  then  it  is  switched  to  *  M£.  The  constraint  (C-3)  gives  a 
measure  of  the  amount  of  trajectory  maneuvering  allowed,  or 
required,  in  order  that  the  vehicl  i’s  position  and  velocity  fall  within 
some  specified  region  at  time  t  =  T. 

Pointwise  constraints  of  the  form  (C-l), (C-2)  specify  the  control 
functions  to  belong  to  a  set  of  U  for  each  tel,  whereas  global  con¬ 
straints  of  the  form  (C-3),(C-4)  specify  the  controls  as  subsets  of 
certain  function  spaces.  In  general,  the  control  functions  are 
considered  to  be  bounded,  measurable  functions  and  the  control 
set  U.  (considered  as  a  subset  of  a  function  space)  is  as  follows: 

0.  =  {u:  u  is  bdd, mble,  u(t)*U(t)  for  t «  i} 

where 

t7  =  {u(t):(C-l)  or  (C-2)  hold} 

or  else 

lL-{u:  u  is  bdd,  mble,  (C-3)  or  (C-4)  hold} 
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For  analytical  purposes,  (C-3),  with  no  pointwise  constraints 
would  be  the  easiest  to  handle.  However,  the  cost  function  should 
then  be  such  that  meaningful  expressions  for  u  are  obtained  when 
the  Hamiltonian  is  minimized.  For  example,  if  it  is  only  required 
to  minimize  Trace  Paar{t)  with  the  constraint  (C-3),  then  the  control 
u  is  eliminated  completely  when  the  Hamiltonian  is  minimized. 

5.4  Cost  Functionals 

Choosing  the  appropriate  cost  function  is  an  important  aspect  of  the 
design  problem.  When  it  is  only  required  to  obtain  good  estimates 
of  die  error  parameters  in  some  fixed  interval  of  time,  then  a 
reasonable  cost  functional  would  be 

J^M)  =  Trace  Wq;F^(T;m)|,  (J-l) 

where  Wa  is  a  positive  definite  weighting  matrix.  More  generally, 
all  of  the  members  of  P(T)  could  be  considered  in  the  minimization 
by  considering  the  cost  functional 

J2m  *  Trace  W P (T  ;  u)]  ,  (J-2J 

t  J 

where  W  is  a  weighting  matrix.  Cu  the  other  hand,  it  may  only 

CSX. 

be  required  that  the  diagonal  elements  of  P  fall  within  a  given 
region  in  the  least  possible  time.  This  could  be  written  as 
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(J-3) 


f  T* 

J3(.«)  ~  I  dt  and  P(T*;m)*G 
•*  Q  ~ 

where  P(T*;m)  denotes  the  elements  of  the  covariance  matrix  P 
which  are  to  lie  within  the  region  G.  The  region  G  is  assumed  to 
have  a  smooth  boundary.  Wim  this  cost  function,  the  question  of 
controllability  must  be  considered.  For  the  minimization  to  make 
sense,  it  must  be  verified  that  controls  exist  which  can  drive  the 
states  P  into  the  region  G. 


It  is  possible  that  it  is  disadvantageous  to  perform  too  much  maneu¬ 


vering  for  the  purpose  of  identifying  the  error  parameters.  In 
this  case  a  term  of  the  form  [**]&(£}{  dt,  which  is  proportional  to 


0 

the  fuel  used,  or  else  a  term  /#?( t)  dt,  which  is  proportional  to 

Jq 

the  energy  used  for  the  maneuvering,  could  be  added  to  the  above 
cost  functionals,,  Thus  we  could  have 


J40/.)  *  Trace  jwp  'T;^+  Cp^T|u(t)|p  dt  p=  1, 2  (J-4) 

JsOO=fTdt+CpfTlM(t)lpdt  p=l,2  (J-5) 

Jo  '  Jq 


T he  weighting  factor  C  might  be  chosen  large  enough  so  that 

r 

pointwise  constraints  on  m  are  not  required.  In  this  case  the  lin¬ 
earization  which  is  made  for  u  in  section  5.2  would  still  be  valid. 
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In  the  next  chanter,  a  particular  optimization  problem  will  be 
specified  by  referring  to  the  various  constraints  (C-i)  and  the  cost 
functionals  (J-i). 

We  can  now  make  the  general  definition: 

The  optimization  problem  is  to  find  a  control  £  such  that 
M  is  a  member  of  the  admissible  class  of  UL,  and 

J  (M)  l  J  (#i)  for  all  m  *  Ll 


5.5  The  Simplest  Examples 

In  this  section  we  consider  the  example 

x  =  px+w,  y  =  mx+ v  ;  cov(w)  *q,  cov(w)*r, 

to  get  an  idea  of  the  minimization  procedure.  In  this  case  the 
covariance  equation  satisfies  the  scalar  equation 

P  “  2pp  ~  r-1mV+q,  p(o)  =  pQ  given. 

Problem  I  J j(m)  =  p(T),  |#(t)|*  Mj  ;  T,MX  given. 

Minimizing  the  Hamiltonian  gives  u  =-M1  sign  (xp),  where  x 
satisfies  the  adjoint  equation 

x  =  -afl/ap  =  (-2p  +  2r~*  m2p)x  ;  x(T)  *  1 
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If  q  >  0,  then  p(t)  >  0  and  v-  -  -Mj  sgn  Aj.  It  is  not  difficult  to 
see  that  x(t)  >  0,  and  therefore  5  -Mj . 

Problem  II  JjjM  =J^  dt,  p(T)  s  a,  and  |ji(t)|  s 

Minimizing  H  =  x-p  +  aq-  l  gives  u(t)  =  -Mi  sgn  (xp)„  At  t  =  T*, 
the  minimum  time,  p(T*)  =  af  H(T*)  =  0,  and  therefore 

A(T)  (-2M!  sgn  A(T)  -r*1  m2  «2  +  q)  +  1  =  0 

Assuming  x(T)  <  0,  means  that  2 lAia  -  2r“*  m2  0^+  q  >0, 
which  implies  p(T*)  is  increasing.  This  contradiction  implies 
x(T*)  *  0  in  10, T],  and  thus  M(t)  =  -Mi,  as  above. 

For  m  =  a,  a  constant,  the  solution  to  the  covariance  equation  can 
be  written  as 

Tanh 

B  or  the  first  problem  T  is  given,  and  is  required  to  minimize 
p(T).  In  the  second  problem  p(T)  =  «  is  given,  and  it  is  neces¬ 
sary  to  minimize  T.  From  this  expression  it  is  clear  that  we 
choose  a  =  -Mi . 
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CHAPTER  VI 


Existence  of  Optimal  Controls  rtnd  Necessary  Conditions 

for  Optimal  Control 

In  this  chapter  the  existence  of  the  optimal  controls  and  the  necessary 
conditions  which  the  optimal  controls  must  satisfy  for  the  optimization 
problems  formulated  in  Chapter  V,  and  for  generalizations  of  these 
problems,  are  considered.  It  is  of  practical  importance  to  insure 
the  existence  of  the  optimal  controls,  since  the  maximum  principle 
of  Pontryagin  is  used  to  obtain  the  necessary  conditions.  Matisx 
notation  is  used  to  specify  the  Hamiltonian,  the  adjoint  equation,  and 
the  equation  which  gives  the  optimal  control  in  terms  of  the  covariance 
matrix  and  the  matrix  of  adjoint  variables. 

Although  the  equations  which  must  be  satisfied  along  the  optimal 
trajectories  are  quite  difficult  to  solve,  and  the  resulting  controls 
are  open  loop  in  nature,  the  results  are  still  of  practical  interest. 

This  is  due  to  the  fact  that  in  the  present  application  to  parameter 
identification,  it  is  entirely  acceptable  to  devote  considerable  effort 
to  obtaining  optimal  trajectory  maneuvers  before  any  experimental 
work  or  ’’flight  testing”  is  performed. 
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6. 1  Existence  of  Optimal  Controls  for  Problem  I 


For  Problem  I  we  consider  a  cost  function  of  the  form 

J l(jl)  =  Trace  WP(T ;  a)  (6.1) 

where  W  is  a  non-negative  weighting  matrix,  T  is  the  specified 
time  of  operation  for  the  identification  procedure,  and  u  is  a 
control  vector.  The  manner  in  which  u  enters  the  covari¬ 
ance  equation 

P  =  A(a)  P  +  PA*(u)  -  PM*  R' 1 MP  +  Q  (6. 2) 

is  assumed  to  be  such  that  the  elements  of  A  which  involve  a 
can  only  be  linear  combinations  of  the  components  of 
u  =  (u\  *  *  • ,  um)*.  Further,  it  is  assumed  that  for  each 
t  e  I  H[to,TJ,u(t)  ^  U(t),  a  compact,  convex  set  in  Em ,  and 
it  is  assumed  that  u  satisfies  the  constraint  (6.3),  where  c  is 
a  constant  vector,, 

u(t)  dt  =  c  (6.3) 

The  control  set  V  is  thus  defined  by  the  set  of  functions 
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convex  set  in  Em,  and  the  integral  constraint  (6.3),  there 
exists  an  optimal  control  u  e  Z/(  such  that  the  cost  functional 
(6. 1)  achieves  a  global  minimum. 


Demonstration  It  is  first  shown  that  there  is  a  uniform  bound 
for  the  elements  of  P,  considered  as  a  vector  P.  That  is,  it 
shall  be  proven  that 


IIPII2  <  ><  » 

(6.4a) 

where  the  norm  1!  •  jlj  is  defined  by 

=  Vtrac*pp  » 

(6.4b) 

and  the  p-  are  the  elements  of  the  matrix  P. 

Indeed,  since 

B  s  and  Q  are  non -negative, 

P  t  A(ii)  P  +  PA*fci)  +  Q 

(6.4c) 

Let  P*  be  a  solution  of  the  equation 

P*  =  APf+  P*A*  +  Q,  P’(to)  =  P0 

(6.4d) 

This  solution  can  be  written  as  (see  Property  7. 1,  Chapter  VII) 
P’(t)  =*(t,t0)P04*(t,t0)+J^t(t,3)Q(s)4.*(t,s)ds,  (6. 4e) 
where  *  (t,  s)  is  a  fundamental  solution  of  the  system 


=  A(t)*(t,s),  =  -*{t,s)A(s);  4(t,t}  =1 

(6.4f) 

By  hypothesis,  the  elements  of  A  will  be  integrable  functions  on 
the  interval  I,  so  that  the  solutions  to  (6.4f)  will  be  unique. 
Further,  it  can  be  concluded  that  P’(t)  is  symmetric  and 


non-negative  for  tel,  since  PQ  and  Q  possess  this  property. 
By  taking  into  consideration  the  problem  from  which  P  (t) 
arizes,  it  can  be  assumed  that  P(t)  is  non-negative.  Thus 

0<<x,P(t)x><  <x,P’(t)x><«  (6.4g) 


For  non-negative  symmetric  matrices 


sup 


<x,Px> 


SUP  !|Jp  x'!2  =  h 
!lx!!=l  1 


where  is  the  greatest  eigenvalue  of  P.  However, 


which  implies  that 


-JU!IPI|2*  HPHi^  l|P|!2  * 

Since  Jj(u)  is  obviously  finitely  bounded  from  below,  we  can 
select  a  sequence  {JikT}€2/  such  that  J(u0  decreases  mono- 
tonically  to  inf  J(u),  where  ueZ/.  We  let  {P^denoie  the 
solutions  to  the  Riccati  equation  (6,2)  which  correspond  to 
the  control  sequence  {  u^,}  .  For  convenience,  Equation (6. 2) 
is  written  in  the  vector  form 

P  =  F(P,u,t)  (6.5a) 


where  the  definition  of  the  vector  function  F  is  derived  from 
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Equation  (6,2).  Since  the  set  U  is  compact  and  convex,  %( 
will  be  weakly  sequentially  compact  for  the  compact  interval 
I  =  [{q,  T] .  Thus  a  subsequence  {u^}  can  be  selected  such 
that  JSk(i)  J“(t)  -or  I«  The  solutions  to  (6.5a)  for  each 
control  Ui  are  given  by 


Pi(t)  =  p0+JT  IteiO>.si<»»»T>dT 


(6. 5b) 


Since  |Fl<  <*iPl  +  $  ,  and  iPi  is  uniformly  bounded,  the  sequence 
{  Pfc}  *orms  a  uniformly  bounded  and  equicontinuous  family  of 
functions.  The  theorem  of  Ascoli  then  assures  us  that  a  sub¬ 
sequence  (Pjj(t)}  converges  uniformly  to  some  function  P(t), 
where 

K t)  -  P„  +£ to  F(PkM,Uk«.-) dT  (6.5c) 


Similarly,  since  Jj  is  continuous  in  P^, 

lim  JiCujj)  =  lim  trace  WPk(T  ;  Ufc)=TraceWP(T) 
k—*»  k  "•  * 

(6. 5d) 


It  is  still  required  to  shew  that  P  is  the  response  to  u  and  that 
u  e  £/.  This  is  accomplished  by  the  techniques  which  are 
used  in  existence  theory.  In  particular,  by  the  assumptions  on 
the  way  u  enters  the  covariance  equation,  F  can  be  written  as 


F(P,u,t)  =  G(P,t)  +  H(P,t)u 


L25 


(6. 5e) 


gag  uiwftffwww 


and  weak  limits  are  taken  to  show  that  P  is  the  response  to  u, 
as  in  page  40  [6 j.  In  this  case  Lebesgue’s  dominated  conver¬ 
gence  theorem  and  a  theorem  on  almost  uniform  convergence 
and  the  continuity  of  G,  H  and  of  the  partial  derivatives  of  G,H 
with  respect  to  the  p^  are  used  to  show  the  weak  convergence 
of  H(Pk(t),t)  to  H(P(t),t)  where  P^  is  the  response  to  u^  and 

A  „  A  — 

P  is  the  response  to  u,  and  P  =  P.  The  compactness  and  the 
convexity  of  the  restraint  set  U(t)  are  used  to  show  that 
u(t)  c  U{t)  for  a.  e.  t  e  I.  Then  u  (t)  is  redefined  on  this  null 
set  so  that  u(t)  €  U(t)  for  all  t  e  I. 

This  implies  that  J(u)=  inf  J  (u)  and  u  is  the  required  opti- 

-  u  tV  ~  - 

mal  control. 

6.2  Existence  of  Time  Optimal  Controls 

The  time  optimal  problem  is  formulated  as  follows.  For  the 
Riccati  equation  (6. 2),  find  a  control  u,  where  u  belongs  to 
some  admissible  set  2i,  such  that  certain  elements  of  P  are 
less  than  some  prescribed  values  in  minimum  time.  In  partic¬ 
ular,  we  shall  require  that  the  diagonal  elements  of  P,  pjj,  be 
less  than  car  j .  To  be  specific,  K  is  defined  by  (6.6),  as  in  the 
previous  section. 
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2/=|u:u(t)  c  U(t),  ^  u(t)dt  =  c|  (6.6) 

Property  6.2  IS.  U(t)  is  a  compact,  convex  set,  then  there 
exists  a  control  ue V  such  that  pjj  s  a j ,  j  =  1,  *  •  • ,  n.  in  mini¬ 
mum  time,  where  the  p^  are  the  diagonal  elements  of  the 
covariance  matrix  P  which  satisfies  the  differential  equation 
(6.2). 

Demonstration  First  it  is  assured  that  there  exist  controls  in 
which  transfer  the  matrix  P  from  P0  to  PT  in  finite  time 
intervals  [to,T]  by  appropriately  choosing  the  positive  numbers 
aj  .  This  assumption  can  be  made  aid  by  solving  (6.2)  for 
arbitrary  ue  2/  and  then  observing  the  p^’s  until  such  time  that 
they  reach  satisfactory  levels  or  j .  The  target  set  for  P-j-,  call 

it  X(T),  will  be  compact  since  the  vector  P  was  shown  to  be 
uniformly  bounded  in  the  previous  demonstration.  We  let  U' 
denote  the  set  of  all  controls  in  'U  which  transfer  P0  to  P-j.  in  the 
time  T"to»  and  let  T-t0  be  the  infimum.  Thus  the  existence  of 
a  control  vie. 'll'  which  corresponds  to  the  time  T-t0  must  be 
demonstrated. 

The  proof  of  this  result  is  similar  to  that  required  in  Property 
6, 1.  Let  {uk*}  be  a  minimizing  sequence  of  controls  from  t/l 
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which  are  defined  on  intervals  [  tQjT^*]  in  such  a  way  that 
the  sequence  {Tkf}  approaches  T  monotonically  from  above. 
The  controls  defined  in  this  way  define  a  sphere  in  L2  ft0,  T  ]  . 
Thus  a  waaMy  convergent  subsequence  {u^}  exists  such  that 
u^  (t)  — ►  u(t)  weakly  in  Ltf.  T] ) .  Thus  we  must  show  that 
ue2/',  which  means  that  u(t)  must  be  in  U(t)  for  teft^T  j ,  u 
transfers  P0  to  X(T)  in  a  time  T-  tQ.  As  in  the  proof  of  the 
preceeding  property  it  can  be  shown  (page  40  16}  )  that  fUU(t). 
the  equibounded  and  equicontinuous  responses  Pk(t)  to  the  con¬ 
trols  u^(t)  contain  a  subsequence  which  converges  uniformly  to 
P(t),  which  is  the  response  to  the  limiting  control  u(t).  It  is 
then  shown  that  P(T)eX(T),  and  that  JJJ  J(u)=T-t0,  so 
that  u  is  the  required  optimal  control. 

6.3  Some  Basic  Properties  of  the  Riccati  Equation 

In  order  to  discuss  the  existence  and  uniqueness  properties  of 
the  Riccati  equation  (6.7a),  it  is  convenient  to  write  it  in  the 
form  (6.7b)  or  (6.7c),  where  the  definition  of  F(P,u,t)  and 
Fij(P,u,t)  is  clear  from  (6.7a). 

P  =  AMP  +  PACu^-PM+R-^MP+Q,  P(to)=P0  (6.7a) 
P  =  F  (P ,  u ,  t)  (6.7b) 

X^FgfP'U't),  i,j  »!,•••, n  (6.7c) 
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Since  the  matrix  A  involves  the  controls  u,  it  can  usually  only  be 
assumed  that  the  Fjj  are  integrable  in  the  time  interval  ( t0,T] . 
The  basic  existence  and  uniqueness  results  for  the  above  systems 
may  be  found  in  [6]  ,  f  7]  ,  the  Appendix  of  (8]  ,  and  (91  .  The 
following  result  can  be  applied  directly  to  the  above  problem. 

Theorem  Assume  that  the  functions  Fy  in  (6.7c)  and  the 

partial  derivative  3Fjj/apk|  are  continuous  in  P,u,  t,  that  is, 
in  the  space  En2+  m+  Then  given  an  initial  point  PQ  =P(t0), 
where  t0  <r  Ic  Ej,  and  a  measurable  control  u,  with  u(t)e  U(t), 
a  set  in  Em,  for  t  c  I,  there  exists  a  unique  absolutely  continuous 
solution  of  (6. 7c)  on  some  subinterval  I*  of  I,  such  that  P(t  }-  Pc , 

K  there  exist  integrable  functions  M(t)  and  K(t)  on  [t^T] 
such  that 

|Fjj(P,u,t)!s  M(t)  and  laFjj  (P^U/Sp^ |s  K(t) 

(6. 7d) 

for  i,  j,k,l  =  1,  •'  **,n,  and  the  solution  F(t)  with  P(tQ)  =PQ  is 
uniformly  bounded,  that  is,  l!P(t)’i<  v< <*  for  tcit^Ti,  then 
this  is  sufficient  to  insure  that  the  absolutely  continuous  solu¬ 
tion  P{t)  of  (6.7)  is  unique  for  the  interval  i\0 ,  T  ] . 


End  of  Theorem 


The  continuity  of  the  functions  Fy  and  9F^/ a  p  in  (6. 7c) 

with  respect  to  the  variables  P,  u,  and  t  is  clear.  They  are 
continuous  with  respect  to  P  Decause  the  right  side  o£  (6. 7c) 
is  quadratic  in  P,  they  are  continuous  in  u  because  u  enters 
the  elements  of  A  linearly.  They  are  continuous  in  t  because 
M,  R,  and  Q  are  assumed  continuous  in  t .  In  order  to  demon¬ 
strate  the  condition  (6. 7d)  it  is  clear  that  it  will  be  sufficient 
to  show  that  the  elements  of  P(t)  are  uniformly  bounded  for 
teft0,T},  that  is,  fp^{t)l?  Y<  *>  for  i,j  =  1,  and 

te  [  t0 , T ] .  In  the  case  that  U('c)  is  a  compact  set,  a  bound  of 
the  form  used  in  the  proof  of  Property  S.  1  can  be  used. 

Sharper  bounds  can  be  obtained  by  considering  the  equation 
(6  7d)  below  (101  instead  of  (6.4c)  where  S  is  a  symmetric 
matrix  which  can  be  chosen  to  lower  the  bound  Pr(t)  of  P(r.). 
The  solution  of  (6.7d)  is  given  by  (6.4e). 

P'  $  (A -SB )  P '  +  P  ’  (A*  -BS  )+(Q  +SBS )  (6.7d) 

Kalman  f  5 1  has  derived  upper  and  lower  bounds  for  P(t)  using 
observability  and  controllability  properties  of  the  linear  esti¬ 
mation  problem,  T^eae  properties,  which  are  defined  next, 
are  for  the  linear  estimation  problem,  and  are  not  to  be  con¬ 
fused  with  similar  notions  which  might  be  introduced  for  the 
control  problem. 
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Definitions:  The  system 

y  =  Ay+  w  ;  cov(w)  =  Q 
z  =My+  v  ;  cov(v)  =  R 
is  said  to  be  completely  observable  on  [t0 ,  T  ]  s 

M(t0,T)=^^(t,T)M,*,(t)R"1(t)M(t)$(t,T)dt 


(6,8a) 
I  if 

(6.8b) 


is  positive  definite  on  I.  The  system  is  said  to  be  completely 
controllable  on  I  i* 


W(to,T )sf  *(t,T)Q(t)$*(t,T)dt  (6.8c) 

•'I 

is  positive  definite  on  I.  The  system  is  said  to  be  uniformly 
completely  observable  (u,  c,  o. )  if  there  exists  fixed  positive 
constants  oj,  cvj,  /3j,  such  that 

0<  arjl  <  Mft-OpOl/JjI  Vt  (6.8d) 

The  system  is  uniformly  completely  controllable  (u,  c.  c. )  if 
there  exists  constants  ovj,  a2,  /?2»  suc^ 

0<»2I^W(t-  cr2,t)<02I  Vt  (6.8e) 

For  the  parameter  identification  problem 
x  =  Ax+Cac*+  w*,  cov(wx)  =  Qx  ; 

«  -  wa,  cov(w»)  ®  Qft  ;  (6.9a) 

z  =  x  +  v,  cov  { v )  =  R  ; 
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we  obtain 


<*> 


0  I 


M  =  [I  I  0]  ; 


■/r‘V  /r*V 

4>aW 


<*>*Q  *<*>*"  + 


qV 


<P*Qa~ 

Qa 


(6.9b) 

(6.9c) 


It  is  thus  required  that  the  integrals  of  the  above  two  matrices 
be  positive  definite,  and  a  small  random  disturbance  wa  must 
be  added  to  the  unknown  parameters  a  to  insure  controllability. 
Assuming  that  R,  Qx  and  Q°  are  positive  definite,  the  positive 
definiteness  of  the  integrals  of  (6.9b)  and  (6.9c)  will  depend  on 
the  amount  of  "maneuvering"  in  the  <t>xand  <i>tt  matrices,  as  was 
observed  for  the  special  cases  which  were  considered  in 
Section  5. 1. 


Kalman  (Lemma  16.9  and  Lemma  16.10  15)  )  has  obtained  the 
following  bounds  ior  ?(t),  under  the  assumption  that  the  linear 
estimation  problem  is  uniformly  completely  observable  and 
uniformly  completely  controllable.  The  upper  bound  requires 
that  P0  bo  non-negative  definite,  and  the  lower  bound  requires 
that  P0  be  positive  definite.  For  t  s  i0  +  a  , 


jw  i(t-cr,  t)+M(t-cr,  t)J4S  P(t)  S  M~1(t-cr,  t)+W(  t-<r,  t)  (6.10) 

Property  6.3  The  solution  P(t)  te  [to,T] ,  of  the  Riccati 
equation  (6.7a)  is  symmetric  and  non-negative  definite. 

Demonstration  Since  R  and  Q  are  symmetric,  P(t)  and  its 
transpose  P*(t)  satisfy  the  same  differential  equation  and, 
since  P0  is  symmetric,  P(to)  and  P*(t0)  both  satisfy  the  same 
initial  conditions.  Assuming  that  suitable  bounds  as  in  (6.7d) 
are  obtainable,  the  uniqueness  of  the  solution  to  (6.7a)  implies 
that  P(t)  =  P*(t)  for  tc  (t0,Tl.  The  non -negative  definiteness 
of  P  (t)  is  plausible  if  it  is  noted  that  the  variance  in  the  esti¬ 
mate  of  a  costate  y*  (see  ~  ..eorem  1,  ( 1 1 )  is  given  by 
E  [(y*,  y (t)  -y (*))]  *  <y*,  P  (t)y*>>  o 


6,4  The  Necessary  Conditions  for  the  Optimization  Problem 

In  this  section  the  necessary  conditions  which  the  optimal  con¬ 
trol  must  satisfy  are  specified.  The  matrix  form  of  the 
Hamiltonian  and  the  adjoint  equations  have  been  previously 
stated  in  [111  .  The  proofs  which  require  only  straightforward 
matrix  manipulations  are  included.  To  be  specific,  the  problem 
formulated  in  section  6. 1  (Problem  I)  is  considered. 
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Property  6.4  If  \Q  is  an  mx  1  adjoint  vector  which  corres¬ 
ponds  to  the  integral  constraint  Jj  u(t)dt  =  c,  and  a  is  a  matrix 
of  adjoint  variables  with  elements  X;  *  which  correspond  to  the 
elements  p.^  of  the  covariance  matrix  P,  i,  j.  =  1,  •  •  •  ,n,  then  the 
Hamiltonian  for  this  problem  is  given  by 

H  =  <Aq, U>  +  Trace j(AP+ PA*  -PM*R“1MP  +  Q]A*[  (6.11) 


Demonstration  As  in  section  6.2,  the  matrix  Riccati  equation 
can  be  expanded  as 


Pij(t)  =F.j(P,u,t),  i, j  =  l,**\n 


(6. 12a) 


and  we  define  a  further  state  vector  by 


x0  (t)  =u(t) 


(6. 12b) 


Applying  the  Maximum  Principle  of  Pontryagin,  there  results 


H  =<A0  ?U>+£  E  P« 

j-l  13  13 


„  jiV*  *  Pin  ‘ 

+  Trace  :  :  ?  : 


=  <Ac,ja>+ Trace  j  : 


Pm’  ’  ’PnnJL^iu’  ’  *  Nm 


~  <^0,u>  + Trace 


<3T  P)A* 


and  (6. 11)  follows  when  the  right  side  of  (6.2)  is  substituted  for 
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Remark  6. 1  If  the  assumptions  of  section  5.  3  are  made, 
then  (6. 11)  may  be  written  as  equation(6. 13c)  below,  where 


(6.13a) 


(6. 13b) 


H  =  x0^+Trace[(PxxAxx4<Pxo'  AaX  ]  +Tr  ace  [  <p°  V  W“*Aax  j 

(6.13c) 

Property  6. 5  Along  the  optimal  trajectory,  the  adjoint  vari¬ 
ables  satisfy  the  differential  equations 


V  (6.14a) 

A  =  -A* A  -  AA  +  M*R‘1MPa  +  aP(M*R'1M)  (6.14b) 


with  the  end  condition 

A(T)  = ‘apcr)tracefWP(T))  =W  (6.14c) 

Demonstration  Applying  the  maximum  principle  of  Pontryagin, 
along  the  optimal  trajectory  gives 
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Next  we  consider  the  term 


S  »  4  ^  ^  t  ^  J  4  ^  1 

J1  lj  J1  nj 

*  *  *  • 

Trace[PA  A  ]=Trace  P 

•  * 

“jn-lj-"  ajnXn) 

*  +  ...  +  + 

+  •  •  •  + 

+  P«1  /C  a4l\,4  +  •••  +  p  y]a,  X  , 

nl4-'  Jl  nj  inn*-'  jn  nj 

Thus  the  matrix  of  partial  derivatives  becomes 

‘^ajlXlj  *  *  *  ^aJnXlj 

H^LlL2l  =  ;  ;  =  AA. 

9P 

•  *• 

z,a ■-!  i  •  •  •  /  a .  A  . 

L  ^  Jx  nj  ^  Jn  njj  (6.15b) 

Next 

Trace  [PM  r“~MPA  ]  =  Trace  [PM’PA  ] 


/ 

SpljXlj  ’  *  *  £  PIjXn j 

1 

Trace  < 

PM’ 

•  • 

•  • 

•  • 

< 

f 

2XjXlj  *  ’  'iCPnj^nj 

I 

i  J 

137 


Since  the  "state”  variables  p-  are  free  at  t-T,  the  correspond¬ 
ing  adjoint  variables  satisfy 

xij<T>=^TIT  (Trace  [WP(T)]) 

along  the  optimal  trajectory.  This  may  be  written  as 


1<T)  =  fejrofc’a  p«+'"+Ew"J  pi»)]=  w- 


Property  6.6  There  exists  a  solution  of  the  adjoint  equation 

i 

(6, 14b)  which  is  unique  and  symmetric  for  t  e  I  ==[  t0 ,  T  j . 
Furthermore,  if  the  weighting  matrix  W  is  positive  definite, 
then  this  solution  is  positive  definite  for  tel. 


Demonstration  The  uniqueness  result  follows  from  the  fact  that 

A(t)  satisfies  a  linear  homogeneous  equation  in  A  .  Thus  (6.14b) 

o 

can  be  written  as  an  n  x  1  vector  equation 
A  =  G(t)  A,  A(T)  =  W. 

The  elements  of  the  n^xn2  matrix  G(t),  which  involve  the 

optimal  covariance  elements  Pjj(t)  and  the  control  u(t),  are 

integrable  on  [  t0 ,  T  ] .  Therefore  Theorem  5. 1,  Appendix  [  8) 

implies  that  the  solution  A(t),  such  that  A  (T)  =  W,  exists  and 

is  unique.  To  show  that  the  solution  is  symmetric,  we  let 
*  -1 

B  =  M  R  M,  and  take  the  transpose  of  both  sides  of  (6. 14b) . 
This  gives 
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[TiA  ]*=  77A*  =  ~  A*A-A*A*  +  A*PB.+  BPA*, 

and  A* (T)  =  W*  -  W.  Thus  A*(t)  and  A(t)  satisfy  the  same 
differential  equation.  Since  the  solution  to  this  equation  is 
unique,  it  follows  that  A*(t)  =A(t),  for  t  c  I.  The  positive 
definiteness  of  A(t)  is  given  by  Property  7. 1,  Chapter  VH 

Remark  6, 2  For  the  system  specified  in  section  5.3, 
equations  (6. 14b)  and  (8. 14c)  become 

J|  /“  -  r 1  (P^A501 4 PXVX-  AXV“X) ;  AXX(t)=0  (6. 1 6a) 
^  Axa=  -AXXAj  +  r'1(PXXAXa* PXV")  ;  axw(t)  =0  (6.16b) 

^Aa“=  -(AtA““<-A"xA,)  ;  a““(T)  -  Wa  (6.16c) 

Property  6. 7  For  the  cost  functional 

J(u)  =  Trace  WP(T;u),  (6.17a) 

an  optimal  control  jl  satisfies  the  condition 

<XQ , u>  +  Trace  [A(u)PA^+  PA*(u)A*]£ 

~  *  (6.17b) 

s  <A0,u>+  Trace[  A(u) PA*  +  PA  (u)A  J 

for  all  u  c  7/,  where  1/  is  specified  by 

t/  =  {u:  u(t)eU(t),  te  I,  JJj2dfc=£}-  (6.17c) 

In  particular,  if  U(t)  is  a  "cube",  that  is 


U(t)  =  {u(t) :  luj(t)|  ?  Mj  for  te  I,  i  =  1,  (6. 17 d) 


(6. 17e) 


then  the  optimal  controls  are  "bang-bang'*,  that  is, 
u^t)  =  ±  for  t  f  I 
This  result,  of  course,  assumes  that  singular  controls  do  not 
arize. 

Demonstration  Equation  (6. 17b)  is  a  direct  result  of  the  defi¬ 
nition  of  the  Hamiltonian  H  and  the  maximum  principle  since 
the  control  u  enters  the  matrix  A  only.  To  obtain  the  "bang- 
bang"  result  (6. 17d),  it  is  recalled  that  u  enters  the  matrix  A  in 
a  linear  fashion.  Thus  that  portion  of  the  left  sides  of  (6. 17b) 
which  involve  u  can  be  written  as 
m  . 

“i  (»o  +  MP»A  >‘»  («.1W) 

i=  1 

Thus  to  minimize  the  Hamiltonian  in  the  case  U  is  defined  by 
(6. 17d),  we  choose 

ill  =  -Mi  sign  (  ^  +  hi  (P,  A ))  (6, 17g) 

Tnisimpli.es  i»i(t)  =±Mi,  provided  that  +  h|(P, A)  *  0.  If 
the  latter  condition  results  for  a  finite  time  interval,  then  uj(t) 
is  indeterminate. 

Remark  6.3  In  the  case  of  Configurations  I  and  n, 

XX  .x<*‘ 
i  A 

ax  ^cycv 

which  results  in  a  control  law 


APA*+PA*A*  = 


0  Ai 
0  0 


pxx  pXa] 
p«x  paoi 


pxx  p 


pQX  paai 


0  0 
A*  0 
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tJWwiwa*'****" 


r(t)  =  -M  sign[x0-v  Trace  (AjP^+^^A^A^+Trac^AiP^^+P^A* /Ca 

(6.18) 

where  i=  I  for  Configuration  I  and  H  for  Configuration  II. 

Remark  6.4  By  considering  certain  monotone  properties  of 
the  elements  of  A  and  P,  it  is  possible  to  obtain  upper  bounds 
on  the  number  of  switches  the  optimal  control  will  have.  In 
particular,  it  is  clear  that  the  diagonal  elements  of  PGa  are 
non-increasing.  If  certain  assumptions  about  the  noise  covari¬ 
ance  matrices  Q  and  R  can  be  made  as,  for  example,  they  are 
also  non -increasing,  then  it  is  plausible  that  Px<*  and  PM 
have  similar  monotone  properties. 

6.5  Unconstrained  Controls 

If  it  is  assumed  that  the  controls  do  not  have  pointwise  con¬ 
straints  and  the  cost  functional  is  of  the  form  (6. 19c),  with  the 
integral  constraint  JJu(t)  dt  =  c ,  then  the  optimal  control  is 
specified  by  (6.19b),  where  the  components  of  h,  and  hj,  are 
determined  by  (6. 19c).  In  order  to  obtain  such  smooth  controls 
and  at  the  same  time  insure  that  the  control  components  uj(t) 

:  are  not  too  large,  it  might  be  possible  to  adjust  the  weighting 

|  matrix  C(t). 

f 

r 
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(6. 19a) 
(6. 19b) 
(6.19c) 


j2<H)  =  2“ jf  “hWBc  +  Trace  [  WP(T  ;H>] 

u(t)=-C4(t)fx0+h(P,A,t)] 

|  Trace  j^APA  +  PA*a]  } 

In  particular,  for  Configurations  I  and  II  the  optimal  control  is 


.-specified  by 


Hi)  =  -4{a0+  tr[(AiPax+PXQA*i)AXX+(AiPaaAax+PaaA+iA1!a )] }  (6. 19d) 


If  pointwise  constraints  are  also  included  for  these  special  con¬ 
figurations  then  the  optimal  control  is  specified  by  (6. 19e)  where 
Mmin  is  side  (6. 19d).  This  result  follows  from  the 

fact  that  a2H(u)/au2  =  c  >  0  • 

I-Mx  if  (t)  ?  -Mx 
pmm(t)if  »ni n  (*>  *  (6.19e) 

+Mj  if  pmin  (t)  *  Mx 


Property  6.8  For  the  unconstrained  control  problem  formu¬ 
lated  above  with  cost  functional  (6, 19a)  and  with  no  integral 
constraint  and  with  W  positive  definite,  there  exists  an  admis¬ 
sible  optimal  control  u  such  that  J  ^  (u)  achieves  a  global 
minimum. 


Demonstration  Choosing  the  control  u(t)  s  0,  the  corresponding 
response  is  P  (t;o)  and  the  cost  is  J  (o)  =  Trace  f  WP (T ;  o)  J  . 

A  finite  bound  for  iiP(t,o);lcc  =  sup  Ipjj(t)!  can  always  be 
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obtained  as  discussed  in  section  6. 3  or  in  the  proof  of  Property 
6. 1.  It  is  next  asserted  that  all  other  controls  u^  which  yield  a 
cost  such  that 

J  (b)  -  J  <£) 

will  satisfy  !i  P(t ;  upL^sliPft ;  0)11^  for  t  <  I.  This  is  due  to 

tiie  fact  that  at  each  instant  of  time  t*  *  I  there  is  the  same 

amount  of  control  available  when  either  u(t)  =Uj(t)  or  u(t)  =o 

were  used  for  t «  [ t0 , tf  3 ,  Thus  if  the  trace  P  (t  ;  Uj)  > 

trace  P(t;  o)  at  t  =  tT,  the  control 

(  o  tQ  s  t «  t' 
u(t)  =  ]  - 

( Ui(t)  t' «  t  SS  T 

would  yield  a  smaller  cost  than  Uj  since  W  >  0 .  Since  the 
norms  in  En  are  topologically  equivalent,  the  assertion  follows. 
It  is  remarked  that  this  result  is  not  necessarily  true  if  the 
integral  constraint  on  u  is  ir.  force,  or  if  W  &  0 .  Further, 
since  all  the  controls  uj  yield  a  finite  cost,  we  have 

J^i)^^i^)4dtiC^!|yi^fc)l*2  dt>  0<  c<  "  *  (6.20a) 

If  there  is  only  a  finite  number  of  controls  with  the  above 
property,  then  the  existence  of  an  optimal  control  is  immediate. 
Otherwise,  a  subsequence  {  u^}  of  these  controls  can  be  selec¬ 
ted  such  that  J{ur}  tends  monotonically  towards  the  infimum  j 
of  J.  Further,  since  the  unit  sphere  in  L 2  is  weakly 
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sequentially  compact,  a  further  subsequence  can  be  selected 
such  that 

u*  (t)  -2U  u(t)  with  u  (  L2(tt0,T]) 

The  responses  {Pjt}  to  the  controls  {  u^t}  will  be  equibounded, 
and  we  can  show  that  they  are  also  equicontinuous.  Indeed,  for 
two  times  t0  §•  t* ,  t"  s  T,  we  have 

t” 

||  Pr( t‘)  -Pi.(tw)!iM  5  f  J||  G(Pit(t) ,  tJIU+IIH  (Pi.(t) ,  t)ui'0Wco}dt 

Since  the  partial  derivatives  9  Pjj/sp^  are  continuous  for 
i,  j,k,i  =  l,***,n,  P  t  En  andu«Em,  the  above  inequality 
becomes 

5  ci  I  t"-tf|  +  C2  ftiiMii(t)ll1dt 

JV  rt„  ,  (6.20b) 

§  Cl  !  t"-t’i  +  c2  |  t"-t'|  (Jt,  Iiui*(t)|l2dt) 

However,  (6.20a)  implies  that  tor  i*  sufficiently  large. 

cjT|jui,(fc)|||  dt  s  j  +  e  ,  05  €  <  «*  .  (6.20c) 

Inequalities  (6.20b)  and  (6.20c)  give  the  required  result 

l|Pi«(t”)-Pr(t’)!!c0  ?  c2  I  r-l’l  +  cj  I  r-fT’  (6.20d) 

Since  the  sequence  {  Pj.}  is  equibounded  and  equicontinuous, 
Ascoli’s  Theorem  insures  the  existence  of  a  subsequence  {Pj,,} 
such  that 

Pi,»(t)  — -P(t)  uniformly  for  each  tel. 

The  fact  that  the  uj*  enter  the  dynamical  equations  linearly  is 
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then  used  to  take  weak  limits  and  show  that  P  is  the  response  to 

2 

u,  the  weak  limit  of  {uit}.  The  convexity  of  ||u(t)||  will  imply 
that  the  cost  functional  is  weakly  lower  semicontinuous.  Thus 
J  (u)  =  j ,  and  u  is  the  required  optimal  control. 

6.6  An  Example 

At  this  point  it  is  instructive  to  consider  the  system  below  as  a 
simple  illustration  of  the  above  ideas  . 

x  =  ax+  ua+  w* ,  cov  (w*)  =  qx  ; 

®  =  o  +  w“  ,  cov  (w“)  =  q“  ; 

y  -•  x  -t-  v,  cov  ( v )  =  r  ; 

Thus 

ra  u*|  fqx  0] 

A  =  ,  Q=  ,  H  =  r,  M  =  [1,0], 

0  OJ  I  0  qa 

and  the  matrix  Riccati  equation 

P  =  AP  +  PA*-PM*R-1MP  +  Q 

yields  the  state  equations 

pu  =  2(apu  +  up21)  - r_1  p^+qx,  pu(0)  given; 

*>12  =  ap12TUp22-r‘lpllp12’  p12(0)  =  0; 

P22  =■•  -r'ip^+q",  p22(0)  given. 

Assuming  the  integral  constraint  JJ  u(t)  dt=0,  we  define  the 

additional  state  equation 
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x0  =  u>  xo(0)=xo(T)  =0 
The  necessary  conditions  specified  by 

A  =  -A*A  -AA+M*R4MPA+  AP*(M*R4M)* 


yields  the  adjoint  equations 

An  =  2(r  pil-a)X11  +  2r  P12*12 


-1 


*12  =  “UXir(r  Pll'a)X12+r  P12A22 

*22  =  -2UX12 

Va0  =0 


Next  the  constraint  lu|  -  ,  and  the  two  cost  functionals 

Jx(u)  =  cap22(T;u)  +  cxpu(T;u) 

J2(u)  =  cap22(T;u)  +  cu^  u2(t)dt 
are  assumed.  For  Jjand  J2  the  boundary  conditions  for  Ay 
are,  respectively 

An(T)=Cx  Au(T)=0 

A^T)  =  0  and  a12<T)=° 

A22<T>  "  c“  A22(T)  *  =“ 
and  the  initial  values  of  the  A„  are  unknown.  For  Jj,  the 

Hamiltonian  is 

Hj(u)  =  XQu  +  Trace  [  (AP+PA*  -  PM*R  -1MP + Q )  A  *  1 

•i  2 

=  XQu+ X1i(2(apii  +  up21)»r  +  + 

"  2X12(ap12+UP22  -r‘lpllp12)  +  A22('j“-  ^PlV 

which  means  that  the  optimal  control  is  specified  by 
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u(t)  =  -M1sign<V2*11P21  +  2A12P22)  . 

Since  Hi£t)=H(T)  dt,  and  9H/3t  =  0,we  can  set 

H(0)  =  H(T)  to  obtain 

.M1jxoj  +  ca(-r“1pi22(T)+Qa=-M|xo-f2x12(0)p22(0)!  + 

+  xn(0)(2ap11(o)  -r_1p  ^(0)  +  (f )  +  x22(0)  +  q“ . 

This  equation  can  be  used  to  eliminate  one  of  the  unknown 
boundary  values  by  considering  the  various  possibilities  for  the 
signs  of  x0  and  xo+2x12(0)p22(0). 

Numerical  results  show  that  the  optimal  control  is  constant  if 
there  is  no  integral  constraint,  and  switches  once  if  the  integral 
constraint  is  in  force.  For  the  cost  functional  J2,  minimizing 
the  Hamiltonian  implies  that  it  is  necessary  to  minimize 
h(u)  h  (xc+  2Anp21+  2x12p22)  +  cu  u2 
Since  02h(u)/ku2  =  cu  >  0,  we  set 

uminH  "(ao+2ahp21  +  X12p22^cU» 
and  the  optimal  control  is  specified  by 

u(t)  =  -Mx  if  Umin  <  'Mx 

=  umin  ^  1  umin1  "  M1 
u(t)  =  *Mi  if  >MX 

If  u  is  unbounded,  then  u(t)  =  umin(t). 
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CHAPTER  VH 


COMPUTATIONAL  ASPECTS  FOR  COMPUTING  THE 
OPTIMAL  CONTROLS 

In  order  to  compute  the  optimal  controls  using  the  necessary  conditions 
of  the  previous  chapter,  it  is  necessary  to  solve  the  two-point  boundary 
t'aiue  problem  associated  with  the  matrix  Riccati  equation  and  the 
matrix  adjoint  equation.  The  initial  condition  for  the  covariance 
matrix  P  is  specified  by  P0,  and,  generaHy,  the  final  values  for  the 
mati  Lx  of  adjoint  variables  A  is  specified  by  at,  which  depends  on  the 
particular  cost  functional  which  is  being  minimized.  Many  people  have 
been  concerned  with  the  numerical  solution  of  similar  two  point  bound¬ 
ary  value  problems.  The  success  of  a  particular  scheme  depends  on 
having  considerable  insight  about  the  solutions  of  the  particular  problem 
being  solved,  since  a  fairly  accurate  guess  of  the  unknown  adjoint 
variables  A0  is  required.  In  fact,  many  published  results  which  show 
successful  iteration  schemes  for  choosing  the  unknown  boundary  con¬ 
dition  for  reasonably  complicated  problems,  start  with  initial  guesses 
in  which  the  cost  functional  is  quite  close  to  being  a  minimum.  The 
nature  of  the  problem  considered  here  is  such  that  considerable  compu¬ 
tational  effort  can  be  justified  for  the  numerical  solution  of  a  particular 
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identification  problem  due  to  the  practical  interest  in  improved  identifi¬ 
cation.  A  particular  example  in  which  improved  identification  of  IMU 
error  parameters  is  very  profitable  occurs  in  the  flight  testing  of 
inertial  platforms  for  purposes  of  obtaining  the  instrument  errors  under 
actual  operating  conditions.  In  those  cases  it  is  necessary  to  "design 
the  experiment",  that  is,  select  a  nominal  trajectory  for  the  flight  test 
vehicle  so  as  to  obtain  as  much  information  about  the  error  parameters 
as  possible  for  a  particular  flight.  The  material  presented  in  the 
previous  chapters  offers  a  systematic  way  to  maximize  this  information. 
The  same  remarks  can  be  made  in  regards  to  inflight  alignment  and 
calibration  of  inertial  measurement  units,  and  for  the  flight  testing  of 
other  devices  in  which  it  is  desirable  to  identify  error  parameters. 

Due  to  the  complicated  nature  of  the  present  optimization  problem,  we 
are  primarily  interested  in  computation  schemes  in  which  it  is  possible 
to  observe  certain  iterations  of  the  computation  and  then  choose  new 
boundary  values  based  on  these  observations.  The  possibility  of  using 
algorithms  which  are  based  on  a  direct  minimization  of  the  cost  func¬ 
tionals  involved  are  not  considered,  although  such  algorithms  are 
suitable  topics  for  further  investigations.  In  particular,  the  e-method 
which  has  recently  been  studied  for  optimal  control  problems  [is] 


could  be  considered. 


HRPPPfPS^SS^PSSPSS^? 


7.1  Iteration  on  the  Initial  Adjoint  Variables,  A0 

For  this  scheme  a  cost  functional  involving  A-j.  and  PT  is  specified, 

for  example,  as  \ 

i 

i 

I(Aq)  =(?T-PT)*Wp(fT-PT)  +  (2x-At)Wx^t-At)  (7.1)  I 

where  the  tiida  denotes  a  final  value  of  P  or  of  A  resulting  from 
an  initial  guess  on  the  adjoint  variables  A0,  and  Wp  and  Wx  are 
non-negative  weighting  matrices.  U  an  integral  constraint  on  the  * 

controls  is  included,  say  of  the  form  u(t)  dt  =  c  ,  then  we  let 
xd(t)  =/Jl(t)  dt,  and  add  a  term  of  the  form  <(So(t)-c),wu(20(t)-cI> 
to  the  functional  (7.1).  Associated  with  this  cost  functional  are 
the  differential  equations 

P  =  AP  +  PA*  -  PBP  +  Q,  P(to)  =  P0  given;  (7.2) 

A  =  -A* A  -  AA  +  BPA  +  APB,  A(T)  =  AT  may  be  given.  (7.3) 

and  the  control  equation  (which  results  from  minimizing  the  Hamil¬ 
tonian) 

Jft(t)  =  h(P,  A,  t)  (7.4) 

The  idea  in  the  minimization  is  to  choose  a  sequence  {A®}  of  initial 
adjoint  variables  so  that  I(A$)  approaches  zero  after  a  reasonable 
number  of  iterations. 
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A  straightforward  method  for  choosing  the  successive  iterations 


is  to  form  31  (AqJ/sAq  and  then  proceed  in  the  direction  of  steepest 
descent.  Since  I(A0)  =  I^Pt(Aq),At(A0)J  ,  we  have 

2W^-Wat+2wA<Sr^at  <7'5) 


The  partial  derivatives  can  be  approximated  by  A£p/AA0  andAATAAc, 
respectively,  where  ABp  and  AA-p  denotes  the  changes  in  the 
elements  of  Bp  and  A-p  due  to  incremental  changes  in  the  initial 
values  of  A.  The  approximate  derivatives  would  be  obtained  by 
solving  (7.2),  (7.3),  (7.4)  as  initial  value  problems  r.(n+l)/2 
times.  If  all  the  elements  of  Bp  are  unspecified,  then  WPs  0  and 
the  resulting  computations  are  considerably  reduced. 


The  main  idea  in  the  steepest  descent  approach  is  as  follows: 

(1)  Make  an  initial  guess 

(2)  Compute  33/3A0|a^  AI/a^q  =  g 

(3 )  Let  A®  =  a^  +  k  •  ,  where  k  <  0 

(4)  If  I(Aq)  <  I(A«>),  return  to  step  (1)  and  continue  this 
process  until  I(A^)  is  sufficiently  close  to  zero. 

There  are  several  methods  available  for  choosing  the  '’gain”  k„ 
The  value  of  k  should  be  small  enough  so  that  the  minimization 
proceeds  in  the  right  direction  and  yet  k  must  be  sufficiently 


152 


large  so  that  the  minimization  does  not  require  too  many 
iterations.  One  method  for  choosing  k  might  be  to  specify 
a  fixed  percent  change  in  the  cost  for  such  iteration,  so 
that  for  n  fixed,  This  results  in  = 

Another  method  is  to  choose  k  such 
that  I(A0  +  kAj)  is  minimized  with  respect  to  k  where  Aj  is  a 
prespecified  matrix.  This  results  in  a  k  specified  by 
-!|Af  :i2/<Af,  grad  I(Af)>.  If  it  is  practical  to  be  in  the  compu¬ 
tation  "loop",  then  it  is  probably  more  satisfactory  to  select  k 
by  observing  the  resulting  changes  in  1^.  Thus  if  I(Ao+1*)  k 
I (A<|>) ,  then  k  is  too  large,  and  if  i(i+1)  is  nearly  equal  to  l(*), 
then  k  couil  be  increased. 


7.2  Iteration  on  the  Final  Covariance  Matrix,  P(T) 

For  particular  problems,  it  might  be  more  efficient  to  iterate 
on  the  final  covariance  matrix  PT  instead  of  on  the  initial  adjoint 
variables  aq.  The  iteration  procedure  is  essentially  as  in  the 
previous  section,  except  that  the  adjoint  and  covariance  equa¬ 
tions  are  solved  with  the  time  variable  reversed  and  the  iteration 
is  with  .  More  specifically,  the  cost  functional  (7,1)  becomes 

I(Bp)  =  !IP0  -Jo'*wp  +  H^o  ~  Ad'wX  (7.6) 
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and  we  let  t=  t  -  t  in  equations  (7.2),  (7.3)  and  (7.4).  This 
method  has  the  advantage  that  the  backwards  covariance  equation 
might  be  less  ’’unstable"  than  the  forward  adjoint  equation. 


7.3  Selection  of  Initial  Values,  a0  (or  PT) 

Since  the  adjoint  equation  is  highly  unstable  it  may  be  impossible 
to  obtain  convergence  in  the  above  iteration  schemes  unless  good 
initial  choices  of  Aq  (or  PT)  can  be  made.  Firstly,  using  results 
of  W.T.  Reid  [15]  we  show  the 

Property  7 , 1  If  the  cost  functional  is  of  the  form 

trace  [WP(T  ;u)],  then  the  solution  A(t)  of  the  adjoint  equation 
will  be  positive  definite  if  W  is  positive  definite. 

Demonstration:  W.T.  Reid  [15]  considered  the  matrix 

differential  equation 

f  =  H(t)  T  +  TK(t)  (7.7) 

and  showed  that  the  general  solution  is  of  the  form  T(t)  =U(t)CV(t), 
where  U(t)  is  a  fundamental  solution  of  U  =  HU  and  V(t)  is  a  fun¬ 
damental  solution  of  V  =  VK,  and  C  is  an  arbitrary  constant 
matrix.  The  elements  of  H  and  K  are  assumed  to  be  continuous 
on  the  interval  I  =  [tQ,  T]e  If  the  elements  of  these  matrices 
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are  assumed  only  fco  be  integrable,  then  the  solution  to  (7.7 )  is 
a  matrix  whose  elements  are  absolutely  continuous  and  (7.7 )  is 
satisfied  almost  everywhere  in  I. 

Under  the  assumption  that  H  =  K*,  it  is  clear  that  if  U (t>  is  a 
fundamental  solution  of  U  =  KU,  then  U*(t)  is  a  fundamental  solu¬ 
tion  of  V  =  VK*,  and  the  solution  of  (7.2)  under  this  assumption 
will  be  T(t)  =  U(t)  CU*(t).  Now  the  adjoint  equation  can  be 
written  as 

A  =  K*A  +aK,  where  K  e  -A+(?M*R~1M)*  (7.8) 

Using  the  usual  notation  *(t,T)  for  the  fundamental  solution  of 
V  =KV,  the  solution  of  (7.8)  is  A(t)  =  *(t,T)  C**(t,T).  Since 
a(T)  =  W,  and  W  is  positive  definite  by  hypothesis,  the  solution 
is  A(t)  =  *(t,T)  W**(t,T).  This  solution  is  positive  definite 
since  *  *(t,T)  always  exists. 

Remark  This  property  could  be  of  practical  value  in 

the  selection  of  the  initial  values  A0.  If  the  solution  fails  to 
remain  positive  definite  for  a  particular  positive  definite  choice 
of  a0,  then  a  new  choice  can  be  made  without  necessarily  com¬ 
pleting  the  solution  of  the  equations. 
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Since  the  optimal  control  can  be  assumed  bang-bang  in  most 
instances,  it  is  feasible  to  test  the  effect  of  various  control 
functions  upon  the  covariance  functions  in  an  open  loop  manner. 
If,  for  example,  useful  upper  and  lower  bounds  on  the  number 
of  switchings  were  available,  the  search  procedure  would  be 
greatly  simplified.  For  the  example  introduced  in  Section  6. 
the  effect  of  the  number  of  switchings  is  shown  in  figure  7. 1. 

In  this  example  at  least  one  switch  is  required  because  of  the 
integral  constraint  J^u(t)  dt  =  0.  These  results  show  that  a 
control  which  switches  once  at  t  =  0. 5  seconds  seems  to  be  the 
optimal.  Once  several  controls  have  been  tried,  the  one  which 
gives  the  best  results  could  be  used  to  solve  the  adjoint  equation 
backwards  (or  else  give  Bp  directly)  to  obtain  a  good  initial 
guess  A^. 


7.4  Discrete  Form  of  the  Covariance  and  Adjoint  Equations 

For  computational  purposes  it  may  be  more  efficient  to  work  with 
the  discrete  form  of  the  covariance  and  adjoint  equations. 

Actually,  in  some  applications  it  is  necessary  to  work  with  dis¬ 
crete  state  and  observation  equations.  The  discrete  form  of  the 
matrix  Riccati  equation  is  given  by 

pk+ 1 =  pk+ 1 "  pk+  lMk+  l'Mk+  lpk+  2  Mk+ 1 +  Rk+ 1 )  Mk+  lpk+ 1 

(7.9a) 

pk+l  =  *k+l,k  pk  *k+ l,k  +  Qk  <7,9b) 

If  Fk  =  Fk(Pk,%)  is  defined  by  Pk+l  -  Pk>  then,  witli  the  integral 
constraint  on  u,  the  Hamiltonian  becomes 

H  =  <A0,uk>+Trace  [FkAj]  (7.10) 

and  the  adjoint  variables  will  satisfy  the  equation 

Vi  -  Ak  ■  -»/3pk  n.ii) 

along  the  optimal  trajectory.  The  optimal  control  uk  is  obtained 
by  minimizing  H.  When  it  is  possible  to  evaluate  the  matr  jz 
DH/d  Pk,  the  adjoint  equation  will  now  be  in  discrete  form.  It  is 
noted  that  in  minimizing  Hk,  a  neat  result  for  the  optimal  uk  is 


not  obtained  as  in  the  continuous  case  because  of  the  inverse 
matrix  operations  involved  in  equation  (7.9).  A  good  approxi¬ 
mation  which  might  be  considered  is  to  minimize  the  continuous 
version  of  the  Hamiltonian  for  the  selection  of  u^. 


For  the  example  presented  in  section  6.6,  the  covariance  equation 
satisfies 

fr  7r7r'  -  C  m  it’  -  d 

rn*  -  d  qa  -  t?2 


k+1 


(7. 12a) 


-k+  i  -k 


A  I, 


and  the  adjoint  equation  satisfies 

[ -l+(ra^)2  -2rna27r’2  (nan')2 

2abr277'2  (2abr(n,+r-n):r'2-2)  -2am,2(*+r-bri) 

(rb*1)2  2rbir,2(7f+r-bn)  -bnrr,2(2(irr+r)-bn) 

»»  m 

In  these  equations  the  variables  are  defined  as  follows: 

a  -  ea^  b  =  -^-[e^-ljujj  ; 

c  =  J*,(k)  d  =  lfe(k)  e  =  P22(k)  ; 

T  =  a(ca+db)  +  b(da+cb)  +  qx  ; 

71  =  ad+  be  ***  3/  *  +  r  . 

For  this  example, with  cost  functional  Jj(u)  of  section  6. 6,  the 
iteration  scheme  is  thus  to  choose  the  sequence  {  jM) }  = 

(X*,  x^( o ), 0 ) , A. 22( o)}*  such  that  !(*(“))  tends  to  zero,  where 


Ak 

(7. 12b) 


(7. 12c) 


158 


I(i'g)  =  ( A  ^(T )  -  c*)2+(  ))2+(X2y  T )  -c0)2 +w*°(X0(T )) 2  (7. 12J) 

For  each  i,  the  vector  of  partial  derivatives  a/^A^Vaxjj^  was 
obtained  by  changing  the  initial  adjoint  vector  by as  shown 
in  the  table  below,  and  then  resolving  equaLons  (7. 12a),  (7, 12b) 
and  computing  the  optimal  control  from  (7, 12e),  four  times. 

u^t)  =  -M  sign  (Xq+  2Ai\P^1+2xi2P22)  ’  (7* 12e) 

« 

Case  Ax0  AX11^X12  AXgAI/AXjj 

1  X  O  O  O  * 

2  OXOO* 

3  0  0  X  0  * 

4  O  O  O  X  * 

7. 5  Future  Problem  Areas 

The  material  presented  in  this  study  suggests  numerous  areas  in 
which  both  practical  and  theoretical  problems  ran  be  found.  A  few 
of  the  problems  which  the  author  feels  might  be  fruitful  areas  for 
future  research  are  mentioned  in  this  section.  This  collection  of 
problems  is  by  no  means  complete,  and,  depending  on  one’s  particu¬ 
lar  practical  or  theoretical  interests,  additional  problem  areas 
would  be  formulated. 

The  techniques  outlined  in  Chapter  n  might  be  used  to  study  gravity 
errors  and  systems  in  which  radar  is  used  to  obtain  observation 
data.  In  the  case  of  strapped-down  IMU’s,  the  elements  of  the 
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computed  transformation  matrix  might  be  used  to  obtain 
additional  observational  data.  As  mentioned  in  Chapter  n, 
angular  velocity  observations  could  also  be  used  to  improve  the 
identification  of  the  error  parameters.  Additional  work  (both 
experimental  and  analytical)  must  be  directed  towards  obtaining 
realistic  probabilistic  descriptions  of  the  environmental  and 
observational  random  disturbances.  To  supplement  the  experi¬ 
mental  data  presented  in  Chapter  IV  additional  specific  error 
parameter  configurations,  such  as  configurations  which  would 
include  accelerometer  bias,  scale-factor,  and  misalignment 
angles,  should  be  considered  to  obtain  additional  insight  into  the 
trajectory  optimization  problem.  Additional  work  on  the  "no- 
noise”  properties  presented  in  Chapter  V  would  also  be  desirable. 

A  meaningful  area  of  research  in  regards  to  control  theory  appli¬ 
cations  would  be  to  study  the  properties  of  the  two-point  boundary- 
value  problem  associated  with  the  matrix  Riccati  equation  and  the 
adjoint  equation.  Theorems  on  the  number  of  switches  the  optimal 
control  makes,  whether  or  not  singular  controls  result,  and 
whether  or  not  the  optimal  controls  are  unique,  would  be  of  great 
value  in  computing  the  optimal  controls.  Numerical  experimenta¬ 
tion  with  direct  and  indirect  methods  for  computing  the  optimal 
controls  should  be  made.  Indirect  methods  (chat  is,  in  terms  of 
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the  adjoint  equations)  have  been  presented  in  Chapter  VH. 
Although  the  resulting  optimal  controls  are  open -loop,  these  are 
of  practical  interest  because  the  maneuvering  associated  with  the 
parameter  identification  problems  of  flight-testing  inertial 
components,  and  inflight  alignment  and  calibration  of  IMU*s 
can  be  precalculated. 

Additional  theoretical  control  problems  arize  when  the  control 
functions  enter  the  dynamical  equations  in  the  complicated 
fashion  suggested  in  section  5. 2  for  Configuration  HI.  In  this 
case  (and  in  other  cases  in  which  mass -unbalance  terms  are 
included)  the  product  of  the  control  and  its  integral  appear  in 
the  dynamical  equations.  Finally,  there  are  systems  in  which 
the  identification  problem  requires  that  a  linearized  minimum 
variance  estimator  be  used.  In  this  case  a  stochastic  optimi¬ 
zation  problem  must  be  considered,  or  else  Monte  Carlo 
studies  would  need  to  be  made  to  compare  the  covariance  which 
is  assumed  to  be  a  solution  of  the  Riccati  equation  with  the  true 
covariance. 
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APPENDIX  "A" 


Discussion  of  the  Correlated  Noise  Model 

We  consider  stationary,  correlated  gaussian  random  processes,  with 
band-pass  power  spectra  of  the  form  shown  in  Figure  A.l.  Such  pro¬ 
cesses  are  typical  of  many  physical  systems  -  in  particular  the  vib¬ 
rations  of  airplane  wings  and  rocket  vehicles.  In  this  appendix  we 
shall  be  concerned  with  the  selection  of  shaping  filters  for  use  in  the 
Minimum  Variance  Estimator  equations,  the  step-size  for  faithful 
numerical  integration,  ana  the  digital  generation  of  the  random  pro¬ 
cesses.  This  appendix  is  included  because  certain  of  the  results 
are  used  in  previous  chapters,  and  it  provides  practical  information 
which  is  not  easily  discerned  from  the  literature. 


Figure  A.l  Experimental  Figure  A. 2  Approximating 
Spectrum  Spectrum 


A.l  Shaping  Filter 

In  order  to  determine  a  suitable  shaping  filter,  the  experimental 
power  spectrum  shown  in  Figure  A.  1  needs  to  be  approximated 
by  an  analytical  function  of  the  form  shown  in  Figure  A.  2. 
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Power  spectrums  which  might  be  used  to  represent  the  given  power 


spectrum  S  (w)  are 

A 

S(1)  (u>)  =  A2  {  exp  -  +  exp  }  (A.  1) 


and 


(«) 


Ag  Ag 

(to  -  a)^+  (co  +  a)^+»02 


(A. 2) 


We  next  approximate  the  given  power  spectrum  (based  on  assump¬ 
tion  and/or  experimental  data)  by  the  fourth  order  rational  function, 
and  then  synthesize  a  shaping  filter  for  this  rational  spectrum 
using  frequency  domain  techniques. 


It  should  be  remarked  that  the  purpose  of  the  shaping  filters  dis¬ 
cussed  here  are  not  for  generating  random  variables.  They  are 
used  to  give  our  problem  the  required  canonic  structure,  and  thus 
the  approximations  made  here  are  not  expected  to  be  as  critical  as 
in  the  case  of  tho  shaping  filters  used  for  generating  the  correlated 
random  process  (Section  A. 2). 

Given  the  experimental  data,  S~(w),  the  fourth  order  approximating 
spectrum,  S  (w),  would  be  obtained  by  minimizing  (A. 3)  using 
standard  approximating  techniques, 

«  =liS~-Sxil,  (A.  3) 
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where  the  norm  l!  •  |i  might  be  of  the  form 


f  °  (Sx<“)  -  Sx(f  dw 


1  *  p  < 


(A.  3a) 


0  <  w  s  w_ 


lS^(w)  -  Sx{w)|  P=°°  (A. 3b) 


Of  course  the  minimization  is  with  respect  to  the  parameters 
A,a,b,  in  Equation  (A. 2).  We  next  demonstrate  several  pro¬ 
perties  of  the  assumed  shaping  filter. 


Property  1  A  shaping  filter  for  the  spectrum  (A.  2)  is  given 


Hx(s)  ='®Ais+c) 
(s  +  b)S  tsr 


(A.  4) 


Demonstration 


S  (s)  * 
xv  ' 


2A2  -  s2  +  b2  +  a2 
s^  +  2s2  (a2  -  b2)  +  {^  + 


Letting 


there  results 


6  =  Tan’1  2ab  .  and  c2=a2+«2  . 


S  ( =  2A  {-s+  c)  (s  -t  c) 
x  (s2  +  cz  erf)  (s'*  +  c2  e~)Qy 
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Figure  A.  3  Pole-Zero  Pattern  of  S^(w) 

The  transfer  function  of  the  shaping  filter  H^s)  is  obtained  from 
the  left-half  plane  poles  and  zeros  of  Sx(s),  since  we  may  write 

Sx(s)  =  |HX(S)|2  Sjs) 

where  S  (s)  =  1  is  the  power  density  spectrum  of  the  white  noise 

W 

input  to  the  filter,  and  !H  (s)l2  =  Hx(s)  H  (~s).  Thus  there 

X  A 

results 

H  (si  = _ '!FAA?..+„sl, _ 

x  (s  +  ce^0/2)(s  +  ce”30/2) 
and  (A,  4)  follows  by  substituting  for  c  and  e, 

li  Equation  (A.  4)  is  written  in  the  time  domain,  we  obtain 
x  +  2bx  +  (a?  +  b^)  x  =  J2A  w  +  n/2~A  cw 


Thus,  the  representation  of  x  as  given  above  is  not  Markovian, 


because  it  requires  that  the  white  noise  w  be  differentiated.  In 
the  following  property  a  method  for  rearranging  the  transfer 
function  of  Equation  (A.4)  is  given  so  that  a  first-order  Markovian 
process  results. 


Property  2 


where 


By  using  the  arrangement  H£(s)  shown  below 
in  Figure  A.4  a  first  order  Markovian  form 
results.  Furthermore,  this  form  may  be 
written  as 


_d 

dt 


>1 

-b  a 

V 

+  _£L. 

r* 

i 

/2. 

-a  -b 

/2. 

-  a 

w  (A.  5) 


x  5E  2y  j  ,  and  A,  a,  b,  c  are  as  defined  above. 


Demonstration  From  Figure  A.4  we  have 


x  =  Xi  + 


*2 


where 

*1  =  X1  +  -Bj  w 

(A.  5a) 

and 

i2  =  -§'x2+  w 

w- 


Figure  A.4 

Rearrangement  of 
the  Shaping  Filter 


Kar 

XI 

1  H^s) 

1 

1  t  * 

H  Vdi 

1 

1 

! 

VC/ 

+  t 

1 

I 

i 

! 

l  . 
JSSsL 

x2 

_Ll  J2A{s  *  c) _ 1 

i  Hxss) 

(s+W  + 

a  r 

! 

J 

x 


The  random  variable  &  has  the  Markovian  property  which  we  are 
seeking,  and  the  constants  a^,  bj,  c^,  a2,  b2,  and  C£  are  chosen 
such  that  x  has  the  same  power  spectrum  as  x.  Now 


From  these  relations  there  results 

c</b^  =  b+ja  and  c2/^2  =  ^  ~  ia 

After  further  substitution  there  results 

bf  =  b%,  ci  =  c2  bj(b-f  ja),  and  bj  =  1  +  j  ^(l±c/b) 

(A.  5b) 
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r^iSGsFjS 


I 


I 


j 

5 

I 

I 


I 


IP 

•! 

s  i 

I! 


Equations  (A.  5b)  imply  that  (choosing  the  negative  sign) 

H£(s)  =  r — - —  +  _ — £ _  =  JjK  s  +  c 

bl s  +  C1  b2s  +  c2  ib^i?  ?T2b7T  c2 


so  that 


K 


=  ±  1  +  ir  <»-'>! 


This  means  that 


w 


Summarizing,  we  now  have  x  =  x  j  +  x2 ,  where 


*1  =  -  (b  +  j  a)  xj  + 


fr w 

K  _ 

BT  w 

bl  andK  are  constants  specified  above. 


X2  =  -  (b  -  j  a)  X2  +  ^ 


(A.  5c) 


From  the  symmetry  in  equations  (A.  5c)  it  is  seen  that  x ^  sothat 

x  =  Xi  +  x2  =  -(b+jajxj  -vo-ja)xj-rK|~i-+  JLl  w 

A  Lbl  blJ 

=  -2Re  ^(b  +  j  a)  XjJ  +  K  w> 
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and 
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-1 

dt 


w 


V 

= 

"0  1' 

o 

X1 

+  J2A 

[  1  1 

x2. 

-c4  -2b 

x2 

■ 

c-2b 

(A.7) 

Demonstration  Let 


X1  =  x  (A.  7a) 

x1=x2+n/2'Aw  (A. 7b) 

x2  =  -2b  x2  -(a2+  b2)  x«  +  J2  A  (c-2b)w  (A. 7c) 


Substituting  (A.  7a)  into  (A.  7b)  gives 

x2  ~  x  ■  v2  A  w 

which  upon  substitution  into  equation  (A.  7c)  gives 

x  Aw=-2bx+2bN/2  Aw-c2x+./2A(c-  2b)  w 
= -2bx  -  c2x  +  A  cw  (A.7d) 

TMs  equation  has  the  vector  form  of  Equation  (A.7) 


Remark  This  procedure  is  certainly  shorter  than  the  method 
given  in  Property  2,  however,  the  substitutions  are 
not  obtained  in  a  straightforward  manner. 
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A.  2  Digital  Generation  of  a  Correlated  Random  Process  Specified 
by  an  Arbitrary  Band-Pass  Power  Spectrum 

A  method  for  generating  random  variables  digitally  is  to  take  linear 
combinations  of  outputs  from  some  ,,standard,,  uncorrelated  random 
number  generator.  In  this  Section  we  discuss  a  method  for  com¬ 
puting  the  coefficients  in  the  linear  -'ombination  of  uncorrelated 
random  numbers.  The  coefficients  are  chosen  so  that  the  resulting 
random  variable  has  a  power  spectrum  which  approximates  the 
given  power  spectrums  in  a  mean  square  sense. 

A  white  noise  sequence  sampled  at  unity  time  intervals  has  the 
spectral  representation 

Wj  =  J  dZw(u>)  (A. 8) 

where 

E  jdaw(cj)l  =  0  and  E  |dZw(u>)|2  =  dw, 
dZw(-u)  =  dZ„(u >), 

E  dZw(<*>)  dZw(u/)  =  0  if  w  *  w' (independant  increments). 

In  order  to  generate  this  correlated  noise,  we  do  not  need  a 
physically  realizable  scheme, i.e.,  we  can  use  future  values  of 


the  random  sequence,  since  these  would  be  available  in  the 
computer  memory. 


Property  4  Suppose  we  assume  that  Xj  can  be  approximated 

satisfactorily  by  a  polynomial  of  the  form 

xj  =  a_N  Wj~N  +  a-N+l  Wj-N+1  + 

- - +a_1w._1+a0wj-  alWj+1  + 

+....  +  aN_1  Wj+N_!+aj+NWj+N  = 

+N 

=  S^wi-k.  <A-9) 

k  =  -N 


It  is  assumed  that  a^  =  a_k  (corresponding  future  and  past  values 
are  equally  weighed),  then  the  coefficients  which  approximate  the 
given  power  spectrum  in  a  least  squares  sense  are  given  by 


cos  katw  du>,  k  =  0, 1, 


,N. 

(A.  10) 


Applying  Equation  (A.  8)  to  (A. 9)  gives 


(A.  10a) 
dZx(w) 
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Thus 


E  |dZx(w)|2  =  ja  +  2S  ^  ajj  cos  k«j2  E  j  dZw(w)j 


Sx(w)  dw  =  j  &Q  +  2^  ajj  cos  kcoj  d w, 


A 1 

Sx(w)  =  j  a0  +  2^  ^  cos  kcv| 
kTi 


(A.  10b) 


where  S„(u)  is  the  given  power  spectrum  of  the  correlated  noise  x. 


The  quantity  in  the  brackets  of  (A.  10b)  is  real,  so  that  we  can  write 


/  N  \2 

Sx (w)  =  |a0  +  2  2.  ak  003 


(A.  10c) 


We  shall  pick  a  mean-square  error  criterion  for  evaluating  the  co¬ 
efficients  ajj.  However,  in  order  to  obtain  simplified  expressions 
for  these  coefficients,  we  approximate  ySx(u)  instead  of  Sx(a>). 
We  then  obtain  the  mean-square  error  over  the  interval  {-*•,  t]  as 


jSx(u))  -  2^"]  cos  fcj|  d« 
k  =  0 


(Note:  The  prime  on  the  summation  indicates  aQ  is  divided  by  two). 


=  f  S^>  d“  +  4  it  ak2/, 

k  =0 


cos  kudw- 
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"  4^  %  f  coskwda)  +  4/^  ]Vaj.aj  f  coskwcos  j^dc^ 

£f0  J-tT  X  b--n 


k=0  j=0 
k*j 


Now  in  actual  fact,  the  power  spectrum  S  (u>)  will  be  zero  after 
some  sufficiently  large  angular  frequency  w  =  n ,  and  the  sampling 
frequency  would  not  be  unity,  but  1/At,  where  At  is  the  time  (in 
seconds)  between  samples.  Thus,  instead  of  cos  kw,  we  have 
coskAtw.  If  we  choose  the  interval.  1-*/ At,  ^Atl  and  the  sampling 
period  At  such  that  ir/At  >n  and  At  is  compatible  with  integration 
requirements  (next  section)  then 


/. 


*/At  \  i  =  k 

coskAtw  cosjAtw  du>  = 

»/At  l  0  j  -  k 


The  error  e  is  then 

•a 


=  f  S^du+s^2  yj3K(u)  cos  kAtwdw 

k  =  0  k=0 


'Tr/At 


Since  we  are  minimizing  e  with  respect  to  the  a^,  which  are 
differentiable  and  unconstrained,  we  obtain 

rn/At _ 

ak  ~  ^  /  yjsx(u) coskAtw dw  =  0  k  =  0, 

“  ^/ax. 


or 


ak=f -T  •\jsx (w ) cos kAtw dw ,  k  =  0,l,...,N  (A.IOd) 

•'A 
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We  are  sure  that  this  is  a  minimum  since  the  second  derivatives 


with  respect  to  a^  are  >  0.  Since  the  integrand  in  equation  (A.  lOd) 
is  an  even  function,  we  obtain  (A.  10)  as  required. 


Property  5 


The  percent  error  in  the  representation  specified 
by  Property  4  is  given  by 


i  -  jb-  pi-  +  >;  a^n  (a. ii) 

where  P  is  the  power  in  the  given  random  process. 


Demonstration  Since  ^S  (w)  is  an  even  function,  it  has  a 


Fourier  Series  expansion 


yjSylv)  =  y?  cos  &?u>,  wej-r,  rj,  (A.  11a) 


where 


cos  -^r-wdw,  k  =  0, 1, . . .  (A.  lib) 

If  we  compare  equations  (A.  lib)  and(A.  10)  with  */ht  =  r,  we  see 
that  =c*k/2»  Thus  ^Sx(w)  =  +  2  ^  a^coskAtu,  and  the 

power  in  the  random  process  is 


k=  1 


rn/At  /'f2  /*7rAt  xr*  2 

P  =  2/  Sx(w)dw=2 /  Sx(w)dw=/  (an+ 2Y]  airCOSkAtw)  dw= 

JO  JO  Jn/bt  u  fcS 


This  error  is  a  function  of  the  number  of  coefficients  used,  N,  and 
the  error  is  inversely  proportional  to  N.  This  monotonicity  of 
the  truncated  Fourier  series  approximation  is  due  to  the  orthogon¬ 
ality  of  the  cosine  functions,.  As  a  percentage,  Equation  (A.  lie) 
becomes 

^H1-^  |¥+Sa*2]i  (A-nf) 

the  required  result  (A.  11). 


Remark:  Suppose  it  is  required  that  the  coefficients  a^  and  N  are 
chosen  so  that  the  total  power  of  X  be  within  a  of  the  power  in  x, 
that  is, 

|P_Jh<  a  (A.ilg) 

^  P  ' 

This  could  be  a  good  criterion  for  determining  the  number  of  co¬ 
efficients  N. 
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Computational  Results 


The  two  power  spectrums  specified  below  by  Equations  (A.  12a)  and 
(A.  12b)  were  considered. 

S(1>M  -  ■£--  jexp  +  exp  j  (A.  •*) 

and 


S^(«) 


(A.  i3a) 


The  corresponding  autocorrelation  functions  are,  respectively, 


R^)(t)  =  2A26"t2  Cf2//2  cos  a r 


(A. 12b) 


-  ^e-b|T| 


cos  a.r 


(A,13b) 


Computational  results  for  S^^(w)  with  a  =  24,0  radians,  o-  =  10.0, 
and  S^(a)  =  1.0,  and  for  2 ^(w),  with  a =24.0  radians,  b  =  5.0, 
and  S^(a)  =  1.0,  are  shown  below.  In  these  tables 
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TABLE  1  SAMPLING  TIME  At  =  v/fi  =  */75  =  ,0427  sec 


s<2>h 

_ _ - _ 

N 

ak 

a 

ak 

a 

0 

. 455257* 

.387541* 

1 

. 196292 

. 14930 

.113449 

. 18522 

2 

113420 

.07231 

-.066344 

. 14067 

3 

-. 106863 

.00398 

-.095744 

.04768 

4 

-.025387 

.00012 

-.045514 

.02683 

5 

-.002864 

.00007 

-.016735 

.02399 

6 

-.002264 

.00004 

.037336 

.00988 

7 

-.001855 

.00002 

.021766 

.00513 

8 

-.001288 

.00001 

-.004489 

.00487 

*  This  coefficient  has  not  beern  divided  by  two. 


The  rwo  sets  of  curves  on  the  following  pages  show  S^(w)  and 

(2) 

S'  '(w)  and  the  approximating  power  spectrums 


■(««+*  £ 


at,  coskAtw 


k“l 


=  1,2, 


for  N  =  10  and  15,  and  At  =  0.02  seconds  (That  is,  the  curves 
corresponding  to  the  coefficients  of  Table  2). 
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From  these  cur  ,-es  we  see  that  the  exponential  power  spectrum, 
S^\w),  is  reasonably  well  approximated  as  a  function  of  a?  by  N  =  10, 
and  very  well  by  N  =  15.  On  the  other  hand,  the  rational  power 
spectrum,  S^)(u>)  is  not  well  approximated,  even  by  N=  15.  Thus 
N  depends  on  t  and  on  the  shape  of  the  given  power  spectrum. 

TABLE  2  SAMPLING  TIME  At  *  C.02  sec 


. . 

S^(«) 

N 

*k 

O! 

a 

0 

.217369* 

, 175520* 

1 

. 184170 

.27880 

. 139382 

.26105 

2 

. 102510 

. 14710 

.061514 

. 18285 

3 

.  013242 

.14490 

-.002265 

. 18066 

4 

-.047379 

.  11677 

-.029463 

. 16222 

5 

-.067106 

.06033 

-.037714 

.13207 

8 

-.056730 

.01959 

-.044446 

.09016 

7 

-.035391 

.00430 

-.044527 

.04810 

8 

-. 170496 

.00065 

-.029811 

.02924 

9 

-.006412 

.00014 

-.008780 

.02761 

10 

-.002095 

.00008 

.004843 

.02711 

11 

-.001062 

.00007 

.010375 

.02482 

12 

-.091053 

.015116 

.01998 

13 

-.001071 

.018987 

.01233 

14 

-.000988 

.016050 

.00686 

15 

-  000845 

.007177 

.00577 

*  Coefficient  has  not  been  divided  by  two. 
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Figure  A.  5a  Approximation  of  Rational  Spectrum 


Approximation  of  Exponential  Spectrum 


A.  3  Step  Size  for  Digital  Generation  and  Integration  of  the 

mm*"**  i  •  «mi».  *  «i  iwnumww— mt*i  ■unwiw  *■—  r  m>mm  »■  n  mm  w36w*,Wt>y»'>WL<fc  VHMMiUt 

Correlated  Hansom  Process 

We  consider  &  zero-mean  stationary,  gaussian  random  process 
with  power  spectrum  S^}  and  auto  covariance  inunction  R(r)  given  fcgp 

R(t)  ~  j  £{w)  coswrdw 

where  it  is  noted  that  here  the  one-sided  spectrum  5  (&?)  is  used- 


A  good  indication  of  the  step  sizes  can  be  obtained  if  the  number  of 
zero-crossings  per  unit  time  the  number  of  maximum  and  minimum 
per  unit  time  and  the  t,  m.  s»  value  of  the  derivative  are  known. 

For  stationary,  Russian,  random  processes  the  expected  values  of 
these  quantities  can  be  calculated  relatively  easily  in  terms  of  the 
power  spe drums  or  autocovariance  functions.  Since  the  random 
process  is  assumed  zero-mean  and  gaussian,  the  probability  dis¬ 
tribution  is  completely  specified  by  R(0)^.  That  is,  the 
probability  that  x(t)  lies  between  x  avu  x+  d:  is 

•7 mJ^Lmmr  exp-X^2R(0). 

^2  “mo) 

Thus  the  probability  distributions  of  such  processes  are  completely 
specified  by  their  variance. 


185 


■gay  w  fflWPjj 


The  expected  number  of  zeros  per  second  at  x\t}  is  §.ven  by 
(Equatteii(3.2-12}  of  tl2]) 

r  *  ~t/2 

r  /Ilf  ?  i  L  w2S|y)i* 

LhE  [No.  of  zeros/aee  «  =  £  1*% - 

t  j  L  Rw  J  \  i  ,  x 

\J0  s(w)dw 

L  (£.14) 

The  expected  number  of  maxima  and  minima  is  given  by  (Equation 
(3.3-S  of 


r 


r*  i 

EjSo.  of  max  and  min  per  sec  ~^r"(oJ 


<rv)  (mJ1^ 


jfw4S(«)dw 


“|l/2 


*  I  Ao 

u?2S(w)dw 


The  r.m.s.  value  of  x(t)  is  given  by 

fl/S 

=  Ixr 

/O 


3/2 


where 


R(")(0)  =  Jfr.  R(r) 


7=0 


(A.  15) 


(A.  16) 


(A.  17) 


Remark  In  the  above  expressions  it  is  assumed  that  the  necessary 
derivatives  of  R(r)  in  the  neighborhood  of  7=0  exist.  According  to 
Theorems  VTE  and  XI  in  section  20. 5  of  A„E.  Taylor’s  Advanced 
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Calculus,  the  necessary  derivatives  oi  R(t)  exist  in  the  neighbor¬ 
hood  of  T=0  provided  that  the  integrals  in  equations  (A.  14  j,  (A.  15), 
(A.  16)  converge.  This  is  a  useful  result  as  we  do  not  need  to  have 
R(r)  explicitly.  For  example,  in  equations  (A.  13a)  and  (A.  13b) 
above,  R2(t)  c*oes  not  have  the  necessary  derivatives  at  r=o,  if 
we  truncate  Sg(w)  all  of  the  necessary  integrals  in  equations  (A.  14), 
(A.  15)  and  (A.  16)  converge,  thus  we  are  sure  by  this  result  that 
the  necessary  derivades  of  H^r)  corresponding  to  this  truncated 
S£(w)  exist.  It  would  be  fairly  tedious  to  obtain  the  R2(t)  corres¬ 
ponding  to  this  truncated  S«j(w) ,  and  then  to  verify  that  the  necessary 
derivatives  exist. 

Generally,  the  relations  (A.  14),  (A.  15)  and  (A.  16)  cannot  be 
calculated  readily.  However,  for  the  exponential  power  spectrum 
(A.  12a)  we  may  calculate  the  quantities  given  by  (A.  14),  (A.  15) 
and  \  16)  in  terms  of  the  derivatives  of  R  j(t)  evaluated  at  t=0. 

In  particular  we  have 

Property  6  For  the  exponential  spectrum,  the  asymptotic 

values  (a  —  »}  may  be  calculated  from 

of  zeros/secj  = 

S  [no,  of  max  and  min/se  cj=  ^ 


E  No. 


2  (A.  18a) 

I  K(0)  *  V 


(A. I8b) 


and 


°i  ■  [Rjt(n)f  =  [-R"(0)]*  -  y'2A<a2+  (X2)  (A.18c) 

Demonstration  Letting 

2  2 

R(r)/2A  =  e“T  a  co8  aT  =  f^(T)  .  f^(r) 
there  results 

R{0)  -  ^(0)12(0)  =  ! 

R"(0)  =  f1(0)f2(0)(,>+2li(0)('>£2(0)(')+£1(0)(")f2(0>  =  -o2-^ 

R(4)(0)  -  £I(0)£^4)+4lj(0)(i)f2(0)!3)  +  6£i(0)<2>l2(o/2)  + 
+4£|0/3^fg(  'fl>  r  ~  a4+6a2tr2+3tr4 

Substituting  these  results  into  {A.  14),  (A.  15)  and  (A.  16)  yields 
equations  (A.  18a),  (A.  18b), (A.  18c),  as  required. 

Experimental.  Results 

Values  of  expressions  (A. 14),  (A. 15),  and(A.16)for  the  two  power 
spectrums  (A. 12a)  and(A.13a),  and  for  various  values  of  the  con¬ 
stants  a,  b,cr,  and  n  are  shown  in  the  curves  below.  The  asymptotic 
values  for  the  exponential  spectrum  (Figure  A. 6)  agree  exactly  with 
the  values  calculated  from  Equations  (A. 18a)  and  (A. 18b)  (that  is, 
when  the  asymptotic  values  are  reached  in  Figure  A.  6. 
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Figure  A.  7  Expected  Number  of  Zero  Crossings,  I  ,  and  the  Expected  Number  of  Maxima 

a 

and  Minima,  I  ,  for  the  Rational  Spectrum. 


To  give  some  physical  Reeling  for  the  above  discussion  and 
results,  a  white  noise  sequence  having  a  flat  power  spectrum  to 
100  cycles/second  was  generated  (white  noise  was  sampled  100 
times/second).  This  sequence  was  than  passed  through  the 
low-pass  filter  (0  to  ’70  radians/second).  By  actually  counting 
the  number  of  zero  crossings  and  the  number  of  maximum  and 
minimum  on  the  time  histories  the  experimental  results  are  com 
pared  below  with  those  obtained  from  the  integral  equations 
(A.  14)  and  (A*  15). 


xe-Interval 

No.  of  Zeros 

No.  of  Max  &  Min 

0-1 

11 

18 

1-2 

14 

18 

2-3 

14 

20 

3-4 

17 

17 

4-5 

14 

19 

5-6 

15 

22 

6-7 

13 

14 

7-8 

13 

17 

8-9 

14 

17 

9-10 

13 

17 

TOTAL 

138 

179 

These  results  show  good  correspondence  between  values  calcu¬ 
lated  from  Equations  (A.  14)  and  (A,  15).  In  particular,  the 
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average  number  of  experimental  zero  crossings  is  13.8,  and 
the  expected  number  of  zero  crossings  from  Equation  (A.  14)  is 
14.0.  The  average  number  of  experimental  maxima  and  minima 
is  17.9,  and  the  expected  number  of  maxima  and  minima  from 
Equation  (A.  15)  is  18.4. 

It  is  important  to  note  that  the  observed  values  of  zero-crossings 
and  maxima  and  minima  do  not  vary  drastically  from  the 
average  values.  This  may  be  readily  observed  by  comparing 
the  columns  of  values  above  with  the  average  values,  13.  8  and 
17. ‘3,  respectively. 


A.4  Comparison  of  Three  Methods  for  Generating 
Digitally  Stationary  Correlated  Bandom  Processes 

In  this  section  we  compare  the  two  models  for  generating  the  cor¬ 
related  random  processes  specified  above  by  Equations  (A.  6)  and 
(A. 9),  along  with  a  third  set  of  equations  to  be  specified  below. 

It  is  instructive  to  make  such  a  comparison  because  each  of  the  . 
three  shaping  filters  considered  here  has  been  formulated  using 
different  assumptions  and  approximations.  Furthermore,  each 
of  the  shaping  filters  would  have  advantages  over  the  other  two  in 
the  application  of  the  minimum  variance  techniques. 

The  first  filter,  Filter  I  (Equation  A.  9)),  is  useful  because  arbit¬ 
rary  power  spectra  may  be  approximated  without  altering  the  filter 
design.  Essentially  this  filter  uses  a  moving  average  of  uncorre¬ 
lated  random  numbers  to  generate  correlated  noise.  This  method 
has  the  disadvantages  that  a  digital  subroutine  is  required  to 
generate  its  coefficients,  and  that  the  filter  is  not  in  canonic  form 
for  use  in  the  error -analysis  equation  of  the  minimum-variance 
estimator.  However,  it  is  useful  for  generating  observations. 

Thi  second  and  third  filters,  Filter  II  (Equation  (A. 6))  and  Filter 
m,  both  require  that  die  gi  ven  power  spectra  be  first  approximated 
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by  a  fourth  order  rational  spectra.  These  filters  use  feedback  to 
generate  the  desired  correlated  noise  in  terms  of  uncorrelated 
random  numbers.  Filters  n  and  III  generally  would  not  represent 
the  given  power  spectra  as  faithfully  as  Filter  I.  However,  these 
filters  are  ir  the  canonic  form  for  use  in  the  error-analysis 
equation.  The  third  filter  probably  gives  a  more  faithful  represen 
tation  of  the  specific  rational  spectrum  considered  than  does 
Filter  n,  since  the  coefficients  of  Filter  HI  have  been  adjusted  to 
improve  the  approximation  in  going  from  a  continuous  random 
process  to  a  discrete  random  process. 

Filter  n  is  used  for  the  same  purpose  as  Filter  in  and  the  deriva¬ 
tion  cf  Filter  n  is  obtained  in  a  straightforward  fashion.  The 
derivation  of  Filters  I  and  n  is  given  above,  and  the  specification 
of  Filter  in  is  given  in  [14  J  .  Combinations  of  the  three  above- 
mentioned  shaping  filters  could  be  useful  in  studying  the  sensitivity 
of  the  minimum-variance  estimates  of  the  state  vector  due  to 
variations  in  the  state  noise  power  spectra. 

We  give  only  a  brief  description  of  the  three  shaping  filters 
considered.  Additional  details  of  these  filters  may  be  found  above. 
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Filter  I  (Section  A.  2) 


We  are  given  a  stationary,  gaussian  random  process,  x,  with  a 
power  spectrum  S^w).  We  assume  that  x(jAt)  4  Xj  may  be 
approximated  by 


2j  =  <&,  ®j>  ■  ^  akwj  -k 


(A.  19) 


where  wj  is  a  white  noise  sequence  and  is  obtained  from  a  "standard" 
uncorrelated  random  number  generator,  and  the  coefficients  a^  are 
constants  which  can  be  computed  from  the  power  spectrum  Sx(w)  of 
Xj.  Specifically,  the  a jr  are  given  by  Equation  (A.  10)  and  N  is 
chosen  to  satisfy  (A,  lig). 


Filter  n  (Section  A.  1) 


We  make  the  assumption  that  the  givon  power  spectrum  may  be 
suitably  approximated  by  the  fourth  order  rational  power  spectrum 


Sx(w)  = 


GT-aF  +  b2  r 


(A  .20a) 


where  B  is  chosen  so  that  Sx(±a;  =  A2.  A  model  for  this  process 


is  specified  by  Equation  (A.  5). 


If  it  is  assumed  that 


the  noise  w(t)  is  constant  over  the  relatively  short  computation 
interval  A,  (in  Section  A. 3  quantitative  procedures  are  given 
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which  insure  that  this  assumption  is  realistic;  the  distinction 
between  a  white  noise  process  and  a  white  noise  sequence  will  be 
discussed  below\  then  the  solution  of  (A. 20b)  may  be  written  as* 


Simplifying  the  notation,  (A.  20)  may  be  written  as 


s  "  An-i-i  +  ^n  wi-i 


(A.  20) 


where  the  definition  of  xj,  An,  and  bn  is  obvious  from  {A.  20;. 
The  ’’white  noise”  sequence  wj  _  j  has  the  property 


=0  and  Covfw^Wj^  E 


f 


where 


( 1  for  i  =  j 
i]  l  0  for  i  j 


*Wong  and  Zakai,  [13]  ,  have  shown  that  one  generally  can’t 
replace  a  stochastic  differential  equation  with  the  limit  of  solu¬ 
tions  of  an  ordinary  differential  equation,  as  the  solutions  of 
the  latter  usually  do  not  converge  to  the  solution  of  the  former. 
However,  the  solution  (A.  20b)  would  converge  to  that  of  (A. 6) 
due  to  linearity  of  the  equation. 
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Filter  in 


A  discrete  model  of  ihe  process  defined  by  (A.  20a)  is  g^n  in  [24 j 


xi  =  blxi~l  ~  ^2xl-2  +  alwi-l  +  a2wi-2 
where  is  as  above,  and  where 

bj  =  2e"bA  C0£5  aA  b2  =  er2bA 

_  r 

al  =  Vs V^l-bg) v»0  ^I+bg rgj  - ^T+b2  +  b[' 


(A«21a) 


a2  =  v» 


+  b2  -bj  +  JT+boT; 


(A.21b) 


a  =  sampling  interval  v>Q  = 


Next  let,  yi  k  x^i  -  {*2/^2)  so  that  xj  =  biX^x -b2yi.i  + aiwj.j. 
These  equations  may  be  written  in  matrix  form  as 


(*i\  =  Al  ”b2\  Ai-l\  /* 

W  \i  0/Vi  j  VI 


wi-l 


(A  .21c) 


Again,  to  simplify  the  notation,  we  write  (A.2Jc)  as 
^  =  Am2i~l  +  bjji  wM 


(A.  21) 


K  a  wW«e  noise  process,  w(t)  has  coyariance  R{t)  a(t-r), 
then  the  corresponding  white  noise  sequence  is 

wj  =  w(t  =  and  Cov  (wj)  =  R(t  =  JA)/& 
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The  covariances  for  Filter  I  and  n  were  assumed  to  be  for  a  white 
noise  process,  and  that  for  Filter  HI  was  for  a  white  noise  sequence. 
Thus  suitable  normalizing  factors  must  be  introduced  into  each. 


Specific  Parameters  Used  for  the  Numerical  Comparison 


In  order  to  make  a  quantitive  comparison  of  the  three  shaping 
filters  discussed  above,  we  took  a  =  24,  b  =  5  and  A  -  5.  For 
these  values  of  a  and  b,  a  sampling  rate  of  a  -  .  02  seconds  was 
selected.  This  value  of  A  was  chosen  according  to  the  criteria 
gi\en  in  Section  A.  3.  For  these  values  of  a,  b,  A,  and  a,  the 
coefficients  in  Equation  (A.  19)  are  given  in  Table  n,  Section  A.  2, 
and  the  matrices  in  Equations  (A.  20)  and  (A.  21)  become 


An  = 


.802586  .417835 

1-.417835  .802586 


) 


brr  - 


rx  .605 172  -.8187311 


A 


m 


,1.0 


o 


'nr 


Simulation  Results 

By  actually  counting  the  number  of  zero  crossings  and  the  number 
of  maxima  and  minima  on  the  time  histories  of  simulated  filters 
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the  results  of  the  three  filters  can  be  compared  with  those  estab¬ 
lished  in  Section  A. 3.  The  observed  values  for  a  5  second  sample 
are  shown  below 


Numerical  Results  (a  =  24,  b  =  5,  A=,G2) 

Number  of  Zeros 
in  5  Seconds 

Number  of  Maxima  & 
Minima  in  5  Seconds 

White  Noise 

129 

170 

Filter  I 

50 

111 

Filter  K 

51 

87 

Filter  m 

4? 

96 

From  Section  A.  3,  the  expected  number  of  zeros  in  5  seconds  is 
approximately  51,  and  the  expected  number  of  maxima  and  minima 
in  5  seconds  is  approximately  115  when  sampled  50  times/second. 
For  white  noise  sampled  every  .  02  seconds,  these  numbers  are 
23.9/second  and  33. 6/second,  respectively. 

Since  we  are  assuming  gaussian  random  variables,  it  is  of  interest 
to  count  the  number  of  points  outside  of  1 1  sigma: 


White  Filter  Filter  Filter 
Noise  i  n  m 

Percent  of  Total  points 
outside  of  t  i  sigma 

33  24  32  28 

4 

i 

* 


I 
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